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ABSTRACT 


This  work  concerns  the  experimental  evaluation  of  the 
self-tuning  controller  (STC)  on  a  pilot  scale  distillation 
column.  The  evaluation  of  the  STC  is  done  by  comparing  the 
performance  of  the  STC  with  that  of  well  tune  proportional- 
integral -derivative  (PID)  controllers. 

The  ST  and  PID  controllers  are  used  to  control  the 
terminal  compositions  of  the  column  when  it  is  subjected  to 
step  changes  in  the  feed  flow  rate.  The  controllers  are 
tested  in  the  following  configurations:  top  composition 
control,  bottom  composition  control,  and  simultaneous 
compositions  control.  For  the  self-tuning  controller,  dual 
ends  control  is  performed  using  both  the  multi  loop  and 
multivariable  ar r angements .  The  performance  of  the 
controllers  is  quantified  by  considering  the  sum  of  the 
absolute  error  values  over  a  specified  time  period  after  the 
feed  flow  rate  is  changed. 

The  experimental  results  show  that  for  all  the 
different  configurations  tested,  the  self-tuning  controller 
is  equal  to  or  better  than  the  PID  compensator.  The  STC  is 
found  to  be  very  robust  in  that  it  can  handle  deterministic 
disturbances  (feed  flow  changes)  as  good  as  a  well  tuned  PID 
compensator  even  though  feedforward  control  action  is  not 
added.  With  the  addition  of  feedfoward  compensation,  the 
STC  performance  is  superior  to  that  of  the  PID  controller. 
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1 .  INTRODUCTION 


1.1  INTRODUCTION 

Interest  in  improving  terminal  composition  control 
strategies  for  distillation  columns  has  increased  in  view  of 
the  dramatic  increases  in  energy  cost.  With  better  control, 
less  over  purification  is  needed  to  guarantee  product 
specification.  The  result  is  reduced  expenditure  in  both 
heating  and  cooling  medias  hence  lowering  energy 
requi remen ts . 

Over  the  years,  considerable  effort  has  been  made  to 
find  better  control  methods  that  would  imporve  the  control 
performance  over  that  of  the  conventional  feedback 
propor t iona 1  - i ntegra 1  -  derivative  (PID)  controller  [1], 

Some  of  the  approaches  have  been  implemented  on  the  pi  lot 
plant  distillation  column  at  the  University  of  Alberta  with 
varying  degree  of  success  [2-9].  The  difficulties 
experienced  with  many  of  the  strategies  that  have  been 
tested  stem  from  the  fact  that  the  techniques  were  designed 
to  handle  linear  systems.  Unfortunately,  the  distillation 
column  is  clearly  a  nonlinear  system  [9].  When  a  controller 
is  tuned  to  perform  well  for  one  disturbance,  the  control 
behavior  for  another  disturbance,  or  set  point  change  may 
not  be  satisfactory.  To  compensate  for  nonlinearity,  the 
control lers  have  to  be  retuned.  An  ideal  controller  would 
be  an  adaptive  controller  where  the  parameters  of  the 
controller  change  to  compensate  for  changes  in  process 
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condition.  The  self-tuning  regulator  (STR)  is  such  a 
controller  [10].  In  the  self-tunning  regulator,  a  recursive 
identification  scheme  is  used  to  identify  either  the 
parameters  of  a  stochastic  model  for  the  process  or  the 
controller  parameters  directly.  In  either  case,  the  control 
parameters  will  be  changed  to  suit  changes  in  operating 
condi t ion . 

The  self- tuning  regulator  has  been  developed  from  work 

o 

done  by  Astrom  and  Wittenmark  [10-11],  Peterka  [12]  and 
Clarke  et  al.  [13-18].  As  originally  proposed,  the 
controller  was  for  a  single  input-single  output  (SISO) 
system.  Borisson  [19],  Keviscky  et  al.  [20-21]  and 
Morris  et  al .  [22-24]  extended  the  work  to  handle  a  class  of 
multivariable  cases. 


1.2  OBJECTIVE 

The  object  of  this  study  is  to  exper imenata 1 ly  evaluate 
the  performance  of  the  self-tuning  controller  (STC)  as 
modified  by  Morris  et  al.  [24]  on  the  pilot  plant 
distillation  column  at  the  University  of  Alberta.  The  STC 
will  be  studied  in  single  loop,  multi  loop  and  multivariable 
configuration.  Results  for  the  self-tuning  controller  will 
be  compared  to  those  for  well  tuned  proportional -integral - 
derivative  (PID)  controllers.  Both  the  STC  and  PID 
controller  will  be  used  to  control  the  terminal  composition 
of  the  distillation  column  when  the  column  is  subjected  to 
feed  flow  disturbances. 


. 
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1.3  THESIS  ORGANIZATION 

The  thesis  will  be  separated  into  six  chapters.  The 
first  chapter  gives  an  introduction.  Chapter  two  contains  a 
survey  of  the  work  done  to  date  and  a  brief  discussion  of 
the  theory  underlying  the  STC.  Extensions  to  the  originally 

o 

proposed  scheme  by  Astrom  and  Wittenmark  [10]  will  be 
presented  and  the  actual  algorithms  used  in  this  study  will 
be  derived.  Chapter  three  describes  the  experimental 
equipment  used  in  this  study.  Chapter  four  deals  with  the 
selection  of  the  manipulated  variable  for  control  of 
distillate  composition.  An  analysis  of  the  internal  reflux 
ratio  will  be  given  to  show  which  of  the  manipulated 
variables,  reflux  flow  or  distillate  flow  should  be  used. 
Experimental  open  loop  responses  will  be  presented  to  show 
the  sensitivity  of  distillate  composition  to  reflux  flow  and 
distillate  flow  manipulations.  Chapter  five  contains  the 
expermental  results.  The  results  will  be  presented  in  the 
following  manner: 

a.  single  loop  control  for  overhead  composition; 

b.  single  loop  control  for  bottom  composition; 

c.  simultaneous  control  in  multi  loop  configuration; 

d.  simulatneous  control  in  multivariable 

configuration. 

Finally,  Chapter  six  summarizes  the  conclusions  and 
recommendations  borne  out  by  this  study. 


2.  DEVELOPEMENT  OF  THE  SELF-TUNING  CONTROLLER 


2.1  INTRODUCTION 

o 

The  self-tuning  regulator  as  proposed  by  Astrom  [10]  as 
the  name  implies  is  an  adaptive  controller.  The  algorithm 
consists  of  two  stages  at  every  sample  interval.  The  first 
stage  is  a  recursive  least  square  identification  scheme  to 
obtain  the  parameters  of  a  linear  stochastic  model  used  to 
describe  the  process  being  controlled  while  in  the  second 
stage  a  control  signal  is  calculated  using  minimum  variance 
of  the  output  as  a  criterion.  The  parameters  estimated  from 
the  identification  step  are  used  directly  as  if  they  are  the 

o 

true  parameters.  Astrom  [10]  has  shown  that  for  systems 
with  constant  parameters,  if  the  parameters  converge,  the 
controller  will  converge  to  the  true  minimum  variance 
controller.  The  structure  of  the  self-tuning  regulator  is 
shown  in  F i gure  2.1. 

There  are  many  advantages  in  using  a  self- tuning 
regulator.  First,  it  is  an  adaptive  controller.  Since  the 
parameters  of  the  model  are  being  constantly  updated, 
changes  in  the  process  parameters  are  automatically 
accounted  for.  For  a  nonlinear  plant,  the  parameters 
identified  will  give  a  linearized  approximate  model  of  the 
nonlinear  plant.  As  the  operating  conditions  change,  the 
parameters  readjust  to  linearize  the  model  around  the  new 
operating  point.  For  a  plant  with  time  varying  properties, 
again  the  parameters  will  adapt  to  the  new  condition. 
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Figure  2.1  Basic  Structure  of  the  Self-Tuning  Regulator 
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Second,  as  it  is  self- tuning,  limited  Knowledge  of  the 
system  to  be  controlled  is  needed.  There  is  no  need  to  do  a 
lengthy  background  identification  test  to  obtain  model 
parameters.  Although  certain  parmeters  have  to  be  speci¬ 
fied,  these  parameters  are  readily  obtained  in  comparison  to 
determining  the  constants  for  a  well  tuned  PID  controller. 
Three,  since  a  model  is  being  used  to  do  the  control 
calculation,  there  is  an  inherent  predictor.  For  a  system 
with  time  delay,  the  predictor  removes  the  effect  of  time 
delay  on  the  system  characteristic  equation.  The  process  is 
controlled  as  if  there  was  no  time  delay.  Finally,  by  using 
a  stochastic  model  formulation,  random  disturbances  such  as 
measurement  or  process  noise  are  automatically  accommodated. 


2.2  LITERATURE  SURVEY 

2.2.1  Historical  Development  of  Self-Tuning  Regulator 

The  idea  of  combining  estimation  and  control  in  one 
algorithm  is  not  new.  Kalman  [25]  in  1958  formulated  a 
simplified  algorithm  that  attempted  to  do  both  estimation 
and  control  at  the  same  time.  However,  the  theory  was 
insufficient  and  digital  computers  were  not  available  for 
process  control . 

Astrom  was  first  involved  with  the  self-tuning 
regulator  when  he  and  Bohlin  presented  a  general  form  for 
stochastic  models  [26].  They  considered  the  maximum 
likelihood  technique  for  identifying  the  parameters  of 
linear  stochastic  models  structured  according  to  their 
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o 

general  form.  Subsequently,  Astrom  [27]  looked  at  the 
accuracy  of  the  estimates  obtained  and  the  possibility  of 
using  recursive  least  squares  to  identify  the  parameters. 

o 

Peterka  [12]  then  used  the  Astrom-Bohl in  model 
formulation  to  derive  an  algorithm  that  combined  estimation 
of  the  model  parameters  and  calculation  of  the  control  law. 
He  used  recursive  least  squares  to  identify  the  parameters 

o 

of  the  Astrom-Boh 1 i n  model;  with  the  model  estimates 
available,  a  minimum  variance  of  the  output  was  used  as  the 
criterion  for  determining  controller  action.  For  some 
conditions,  Peterka  found  that  the  controller  used  was 
identical  to  the  minimum  variance  controller  that  could  be 
derived  had  the  model  parameters  been  known. 

o 

Following  Peterka' s  work,  Astrom  and  Wittenmark  [28] 
did  a  thorough  study  looking  at  seperate  and  combined 
estimation  and  control  stratergies.  They  considered  models 
that  have  constant  unknown  parameters  and  models  that  are 
time-varying  Gauss-Markov  processes.  With  Peterka' s  results 
and  their  own  work,  Astrom  and  Wittenmark  formulated  the 
self-tuning  regulator  [10].  They  proved  two  theorems 
concerning  the  self-tuning  regulator.  Using  weak 
assumptions,  they  showed  that  if  the  parameters  converged, 
then  the  covariance  of  the  output  variable  and  the 
crosscovar i ance  between  the  output  and  input  variable  would 
tend  to  zero.  Futhermore,  it  was  shown  for  systems  with 
constant  parameters  that  if  the  parameters  converged,  the 
controller  would  converge  to  the  minimum  variance  controller 
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that  could  be  derived  if  the  model  parameters  were  Known. 
Results'  of  a  detailed  simulation  study  using  the  self-tuning 
regulator  were  presented  by  Wittenmark  [29]. 

A  more  complete  survey  of  the  different  studies  that 
led  to  the  formulation  of  the  self-tuning  regulator  has  been 
given  by  Chang  [30].  The  papers  by  Clarke  [18]  and 

o 

Astrom  [11]  are  excellent  reviews  of  many  of  the  self-tuning 
studies  that  have  been  performed. 

2.2.2  Modification  and  Extension 

In  its  original  formulation,  the  self-tuning  regulator 
had  very  few  options.  For  instance,  there  was  no  set  point 
following  capability.  Many  extensions  and  modifications 
have  since  been  added  to  increase  the  generality  and 
robustness  of  the  controller. 

Wittenmark  has  extended  the  regulator  by  adding 
features  to  allow  set  point  following  for  the  output  varia¬ 
ble  and  feedforward  action  for  measurable  disturbances  [29]. 

Clarke  et  al .  [13-18]  have  presented  a  generalized 
version  of  the  self-tuning  regulator.  The  term  self-tuning 
controller  designation  is  used  for  the  Clarke  approach  to 
distinguish  this  type  of  controller  from  the  self-tuning 
regulator  for  the  Astrom  approach.  In  Clarke's  scheme, 
instead  of  identifying  the  parameters  of  a  model  and  then 
subjecting  the  identified  model  to  some  control  criterion, 
the  controller  parameters  are  identified  directly.  In 
addition,  the  cost  function  for  the  controller  is  modified 
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to  include  control  effort,  process  output  and  set  point 
following.  Inclusion  of  control  effort  into  the  performance 
criteria  achieves  two  objectives.  Adding  a  penalty  for 
control  effort  causes  the  control  output  to  remain  within 
the  constraint  limits.  The  'bang-bang'  controller  response 
that  is  observed  with  minimum  variance  controller  is 
prevented.  Second,  it  can  be  shown  that  with  the 
appropriate  choice  of  control  weighting,  systems  that  are 
nonminimum  phase  can  be  made  closed-loop  stable  [14].  In  a 
treatment  of  nonminimum  phase  systems,  Cegrell  and 
Hedqvist  [31]  have  proposed  a  similar  solution  to  that  used 
by  Clarke.  Astrom  [3.2],  Astrom  and  Wittenmark  [33],  Turtle 
and  Phillipson  [34]  and  Peterka  [12]  have  suggested  the  use 
of  suboptimal  stratergies  to  handle  such  systems. 

Morris  [22-24]  inturn  has  generalized  the  work 
presented  by  Clarke.  He  has  introduced  a  discrete 
compensator  in  PID  form  for  the  control  weighting  used  in 
the  cost  function.  The  PID  controller  is  suggested  because 
no  offset  results  when  a  set  point  change  is  made  and 
futhermore,  the  characteristics  of  the  PID  algorithm  are 
well  understood.  Set  point  model  following  is  added  to 
allow  the  process  output  to  change  according  to  any  user 
specified  response  when  a  set  point  change  is  made;  e.g., 
the  output  can  be  made  to  follow  a  first  order  response  with 
some  specified  time  constant.  By  allowing  the  output  to 
change  slowly,  there  is  no  need  for  excessive  control  action 
at  the  start  of  a  set  point  change. 
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Morris  has  further  extended  the  controller  to  handle 
multivariable  systems  with  an  equal  number  of  inputs  and 
outputs.  The  multivariable  system  is  reduced  to  a  series  of 
decoupled  single  loop  controllers.  The  decoupling  is 
achieved  by  using  the  feedforward  option  to  handle  the 
interaction  between  the  individual  loops;  i.e.,  the 
interactions  between  the  loops  are  treated  as  measurable 
disturbances.  This  approach  removes  the  need  to  identify 
directly  the  parameters  (matrices)  of  a  multivariable 
system.  Direct  identification  of  such  system  has  been 
considered  by  Gustavsson  et  al.  [35]. 

o 

Astrom  and  Peterka  have  considered  extending  the 
self-tuning  regulator  to  the  multivariable  case  for  systems 
with  constant  and  unknown  parameters  [36].  They  looked  at 
identifying  the  parameters  of  the  multivariable  model 
directly.  Borrison  [19]  and  Keviczky  [20]  have  presented 
algorithms  to  handle  multivariable  systems  that  identify  the 
parameters  of  the  multivariable  structure  directly.  De 
Keyser  [37]  has  proposed  a  multivariable  self- tuning  scheme 
that  calculates  a  control  output  based  on  the  identified 
estimates  and  the  covariance  of  the  estimates.  With  this 
approach,  excessive  control  is  prevented  when  the  parameters 
are  bad;  in  effect,  a  more  cautious  controller  is  devised. 

Different  identification  schemes  have  been  considered 
for  the  self-tuning  regulator.  Aside  from  normal  least 
squares,  other  methods  that  have  been  tried  include: 
stochastic  approximation  [38],  extended  least  squares  [39] 


and  maximum  likelihood  [40].  Recently,  it  has  been  shown 
that  the  standard  least  square  identification  method  can 
lead  to  problems  when  a  system  experiences  no  disturbances 
over  a  long  period  of  time  [24].  Without  any  disturbance, 
the  identification  scheme  receives  no  new  information  and  as 
a  result,  the  covariance  matrix  tends  to  'blow-up'  or  lose 
positive  definiteness.  In  either  case,  the  parameters  will 
diverge  causing  stability  problems  with  the  controller.  To 
overcome  this  problem,  Morris  [24]  has  proposed  modifying 
the  normal  least  square  scheme  by  simply  updating  the  square 
root  of  the  covariance  matrix.  Morris  also  proposed  to 
employ  a  recursive  learning  estimator  [24]  once  reasonable 
estimates  have  been  obtained  with  normal  least  squares. 
Clarke  [18]  has  reported  alternate  solutions  which  modify 
the  identification  after  a  'no  disturbance'  situation  is 
detected.  One  possibility  would  be  to  stop  the  identifi¬ 
cation  if  the  determinant  of  the  covariance  matrix  exceeds 
some  specified  limit  or  to  introduce  a  disturbance  if  this 
condition  occurs.  The  third  approach  Clarke  suggested  was 
that  the  forgetting  factor  be  adjusted  so  that  it  tends  to 
unity  if  no  disturbance  is  entering  the  system. 

2.2.3  Appl icat ion 

Numerous  applications  of  the  self-tuning  regulator  for 
the  control  of  various  different  types  of  systems  have  been 
reported.  Those  considered  pertinent  to  this  work  will  be 
briefly  discussed.  Wittenmark  and  Borisson  [40]  have  used 
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the  self-tuning  regulator  for  controlling  the  moisture 
content  of  a  paper  machine.  They  found  that  the  regulator 
was  able  to  handle  nonstationary  disturbances  and  that  the 
choice  for  the  number  of  parameters  in  the  model  was 
unimportant.  Cegrell  and  Hedqvist,  on  a  different  paper 
machine,  used  a  self-tuning  regulator  with  five  inputs  and 
one  output  [41].  They  reported  that  the  self-tuning 
regulator,  in  comparison  to  a  PI  regulator,  reduced  the 
output  variance  during  quality  change  and  mill  step-up.  In 

addition,  their  experience  was  consistent  with  that  of 

Borisson  and  Wittenmark  in  that  they  also  concluded  that  the 
number  of  parameters  in  the  model  was  unimportant. 

Borisson  and  Syding  have  reported  the  use  of  a 
self-tuning  regulator  on  an  ore-crusher  [42].  They  found 
that  the  regulator  was  able  to  adapt  to  changes  in  the 

characteristics  of  the  ore  feed  and  the  crusher  itself.  As 

a  direct  result  of  using  the  self-tuner,  a  10%  increased  in 
production  was  achieved  over  that  possible  using  PI  control. 

Dumont  and  Belanger  [43-44]  have  used  the  self- tuning 
regulator  for  controlling  a  titanium  oxide  kiln.  They  were 
able  to  reduce  the  number  of  parameters  required  in  the 
self-tuning  regulator  by  applying  the  regulator  on  a  plant 
that  was  already  controlled  by  a  fixed  parameter  regulator. 
They  reported  that  with  this  type  of  stratergy,  instead  of  a 
normal  10-12  hours  of  of f speci f icat ion  product  during  a 
grade  change,  the  switchover  time  was  reduced  to  2  hours. 

Keviczky  et  al .  [21]  and  Westerlund  et  al.  [45]  have 
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applied  the  self-tuning  regulator  to  control  of  cement  raw 
material  mixing.  The  application  by  Keviczky  was  on  a 
multivariable  system  and  the  control  criteria  was  based  on 
minimum  variance  of  the  output  and  some  required  average 
over  a  time  duration.  In  both  applications,  it  was  found 
that  the  STR  gave  results  that  were  superior  to  those 
achieved  with  a  conventional  control  approach. 

Various  authors  have  presented  results  from 
applications  of  the  self-tuning  regulator  to  pilot  scale 
equipment.  Sastry  et  al.  [46]  reported  on  the  performance 
of  the  self-tuning  regulator  for  controlling  the  overhead 
composition  of  a  binary  distillation  column.  In  general, 
the  self-tuning  regulator  performed  better  than  a  PI 
controller.  Chang  [30]  used  the  self-tuner  on  a  pilot  scale 
evaporator.  He  considered  how  the  various  design  parameters 
affected  the  performance  of  the  regulator.  Harris  et 
al .  [47]  studied  the  control  of  a  pilot  scale  packed  bed 
reactor  with  the  self-tuning  controller.  The  reaction 
involved  was  extremely  temperature  sensitive  which  resulted 
in  highly  nonlinear  and  nonminimum  phase  system  behaviour. 
The  self-tuning  controller  was  reported  to  be  superior  to 
the  conventional  PI  algorithm.  The  study  also  looked  at  the 
problem  of  '  par ameter -wi ndup'  ;  this  is  a  situation  when  the 
parameters  starts  to  diverge  because  the  controller  output 
is  clamped  at  one  of  the  constraints. 

Applications  of  the  self-tuning  regulator  have  also 
been  reported  for  the  control  of  an  enthalpy  exchanger [ 48 ] , 
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di ges ter [ 49 ] ,  absorption  column[50]  and  ship  steer i ng [ 5 1 ] . 

2.2.4  Convergence  and  Stability  Analysis 

Ljung  has  presented  results  for  convergence  analysis  of 
general  recursive  stochastic  algorithms  [52-56].  His 
approach  was  to  associate  the  algorithms  with  an  ordinary 
differential  equation  so  that  the  convergence  properties  of 
the  algorithm  could  be  obtained  by  looking  at  the  stability 

o 

properties  of  the  differential  equation.  Astrom  has  used 
L jung' s  results  to  analyze  the  self-tuning  regulator  [11]. 

For  the  self-tuning  controller,  convergence  analysis  is 
based  on  the  idea  of  treating  the  controller  as  a  feedback 
system  [57,58].  Using  the  Martingale  theorem  [59,60], 
Gawthrop  [61]  has  shown  that  under  certain  conditions, 
convergence  will  occur  with  a  probability  of  1.0. 

2.3  DERIVATION  OF  THE  SELF-TUNING  CONTROLLER 
2.3.1  Control  Law  Formulation 

A  single  input-single  output  (SISO)  stochastic  process 
can  be  described  by: 

Ay^  =  z"kBut  +  +  z"*<1"Lvt  (2.1) 

where  A , B , C , L  are  polynomials  in  z"^  with  a0  =  1 ,  b0  =  0  and 
cQ  =  1  with  z“*Vj.  denoting  the  value  of  y,  N  samples  in  the 
past.  The  terms  y-j-,  u 5^,  v-^  are  respectively  the  process 
output,  control  input,  random  disturbance  and  measurable 
disturbance.  The  system  is  considered  to  be  of  order  n  with 


J-  .  '  M  ■  ^ 

■ 


•  .  i  !  :  v  *  5.  n  3  y  -9 .  » lj  jd'l  :=|i| 


15 


a  time  delay  equal  to  K  times  the  sample  time  and  kL  the 
time  delay  associated  with  the  known  disturbance.  It  is 
also  assumed  that  C  has  roots  that  are  within  the  Z-domain 
unit  circle  and  that  £  is  a  random  uncorrelated  noise 
signal  with  zero  mean.  The  controller  is  designed  to 
minimize  the  cost  function: 

J1  =  E{  <  pY  t  +  k  '  Rwt)2  +  <Q'ut)2}  (2.2) 

where  P , R , Q'  are  polynomials  in  z'1  which  allow  weighting  to 
be  applied  to  the  output,  set  point  and  control  effort.  The 
term  (  P  y  t  ^  ~  Rw  provides  set  point  following  while  (Q' u-jO 
puts  a  penalty  on  control  action.  From  the  form  of  the  cost 
funtional  it  can  be  seen  that  if  Q'  u t  does  not  tend  to  zero 
at  steady  state,  Py^+k  _  Rwt  will  tend  to  some  constant 
other  than  zero  and  the  process  output  will  exhibit  an 
offset  (cf.  equation  (2.23)). 

In  order  to  minimize  ,  an  expression  is  required  for 
y  .  in  terms  of  past  input  and  output  data.  Here,  y  , 

t+k  K 

denotes  the  k-step  ahead  predictor  for  the  output.  The 
basic  derivation  of  the  self~tuning  regulator  was  given 
first  by  Astrom  [10].  However,  the  ensuing  derivation 
follows  closely  the  formulation  given  by  Clarke  [13]  and 
Morris  and  coworkers  [24]. 

Expressing  Equation  (2.1)  in  terms  of  the  k-step  ahead 
predictor  gives: 


B 


y 


t  +  k  A 


c 

+  zk-S.  + 
A  1 


(2.3) 


( 
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In  Equation  (2.3),  at  time  t,  providing  kL  is  larger 
than  k,  only  the  future  random  disturbances  from  time  t+1  to 
t+k  are  unknown.  Separating  zMc/A)^  into  disturbances  up 
to  the  present  time  and  those  at  future  time  gives: 


z 


=  E'zk?t  + 


(2.4) 


where  E' z^£  represents  the  future  disturbances  and  (F'/A)?t 
equals  the  disturbances  up  to  the  present  time.  The 
polynomials  E'  and  F'  have  terms: 


E'  =  1  +  e^z  +  •  •  ■  • 

.  .  .  +  e-  z-(k-l) 
k-l 

(2.5) 

F'  =  f0  +  fiz  +  •  •  • 

•  •  •  +  f '  z-  (  n  - 1  ) 

(2.6) 

n- 1 

Substitution  of  Equation  (2.4)  into  the  predictor 
equation  (2.3)  and  noting  that  z^S^  is  the  same  as 
gives : 


+  E'C 


t  +  k 


(2.7) 


The  random  disturbance  ^  up  to  time  t  can  be  expressed 
in  terms  of  previous  inputs  and  outputs  by  rewriting 
Equation  (2.1)  in  terms  of  as: 


5 


(2.8) 


Using  Equation  (2.8)  to  substitue  for  5  in  Equation 

V 

(2.7),  after  some  rear r angement  yields: 
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yt  +  k  = 


1 

A 


z"kF/ 


F# 


zk-kLLv^ 

(2.9) 


From  Equation  (2.4),  it  follows  that  the  terms  in  the 
brackets  simplify  to  E'  /C  so  Equation  (2.9)  becomes: 


E'  B  F'  E'  L 

yt  +  k  =  E'  ?t+k  +  -  ut  +  -  yt  +  -  zk-kLv.  (2.10) 

C  C  C  r 


Since  has  been  assumed  to  be  a  white  noize  sequence, 
the  future  disturbance  is  expected  to  be  zero.  Taking  the 
expectation  of  Equation  (2.10),  the  predicted  output  is 
given  as: 


* 


y 


t+k 


E'  B 


E#  L 
C 


(2.11) 


where  the  asterisk  denotes  a  predicted  quantity. 

All  the  terms  on  the  right  hand  side  of  Equation  (2.11) 
are  known  at  time  t  provided  k*-  is  equal  to  or  greater  than 
k;  i.e.,  the  time  delay  in  the  feedforward  path  must  be 
larger  than  the  time  delay  for  the  control  action.  If  the 
above  condition  is  not  true,  a  predictor  for  future  measured 
disturbances  is  needed.  However,  for  most  actual  processes, 
the  time  delay  for  the  measured  disturbance  is  much  larger 
than  the  time  delay  for  the  control  variable. 

The  performance  criteria  given  by  Equation  (2.2) 
requires  a  weighted  output.  Applying  '  P'  weight  to  Equation 
(2.3)  yields : 


■  ...  ■  -  ,  ; .  c  ‘  Hg' 
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PB  PC 

—  u  +  - 

A  r  A 


.  .  L  PL 

Pyt+k  r  r  ut  +  “5t+k  +  zk'k  -  vt 


(2.12) 


Now  by  defining: 


P  =  P  IP 
rN  '  D 


(2.13) 


and , 


PC  z"kF 

-  =  E  +  - 

A  APj 


2.  14) 


Using  Equations  (2.12)  to  (2.14)  and  following  the 
derivation  used  for  the  k-step  predictor  (Equation  (2.10)), 
the  weighted  output  Py^  +  ^  can  be  written  as: 


Py 


t  +  k 


=  EC 


t+k 


F  EB  EL  l 

+  -  v  +  —  u  +  —  zk  ”  k  V 

CP  1  c  t  c  < 


(2.15) 


As  before,  the  best  estimate  of  future  disturbances  is 
zero  so  the  k-step  ahead  weighted  output  predictor  is  given 
by: 


(Py 


t  +  k 


F  EB  EL  .  L 

)*  =  —  y  +  —  u  +  —  zk"kv 

CP  t  c  1  c  1 


(2.16) 


The  actual  weighted  output  and  the  predicted  output  are 
related  by: 


Py 


t+k 


=  (Pyt+k)  +  e 


t  +  k 


(2.17) 


where,  et+k  =  ECt+k-  Substituting  the  weighted  output  given 
by  Equation  (2.17)  into  the  cost  function  yields: 


j,  =  e<  ( (Pyt+k**  +  et+k  -  Rwt>2  +  to'utl2  > 


(2.18) 


. 
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Since  it  has  been  assumed  that  the  random  disturbances 
are  uncorrelated  with  either  the  input  or  output,  the 
expectation  simplifies  to: 

^  [(pyt+k)*  '  Rwt^2  *  +  E{[Q/ut]2}  +  a2+k  (2.19) 

where  is  the  variance  of  the  random  disturbances. 

The  cost  function  is  minimized  by  setting  the  partial 
derivative  3d.j/3u  to  zero.  To  obtain  the  true  optimal 
solution,  the  expectation  operator  must  be  included  in  the 
differentiation  [62],  However,  a  complicated  control  law 
results  so  instead,  a  stage  by  stage  optimal  controller  is 
used.  In  this  work,  following  the  approach  of  Clarke  and 
Gawthrop  [14],  the  expectation  is  removed  from  Equation 
(2.19)  and  the  resulting  equation  is  differentiated  with 
respect  to  u  .  The  controller  that  results  is  optimal  each 
time  the  control  calculation  is  performed.  Although  the 
solution  is  not  globally  optimal,  Harris  [47]  has  shown  that 
the  control  response  of  the  suboptimal  stratergy  is  close  to 
the  result  that  is  obtained  when  the  true  optimal  controller 
is  used.  Removing  the  expectation  operator,  the  partial 
derivative  becomes: 


ad,  *  MPy 

— 1  =  2[ (Py  )  -  Rw  ]  *  - 

8u  aut 


tn)  .  ... 


+  2Q'ut 


(2.20) 


3U. 


Performing  the  necessary  differentiation  and  noting  that 

( Py  )*  is  related  to  u  (cf.  Equation  (2.16)),  setting 
t+k  t 

Equation  (2.20)  to  zero  gives: 


■ 
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l  (Pyt+j<)*  -  Rwt  ]  *  +  qo'Q/ut  =  0  (2.21) 

C 

0 

where  eQ ,  bQ  ,  cQ  and  q0'  are  the  leading  coefficients  for 
the  polynomials  E,  B,  C  and  Q'  .  Defining: 


Q 


(2.22) 


and  recalling  that  e0  and  c0  are  equal  to  1,  Equation  (2.21) 
can  be  writ  ten  as : 


•  ★ 

ut  =  -  (  Rwt  -  (Pyt+k)  ]  (2.23) 

Q 

Examination  of  Equation  (2.23)  shows  that  it  is  similar 
to  a  standard  feedback  control  law  with  the  main  difference 
that  the  feedback  value  is  a  predicted  output  instead  of  the 
actual  measurement.  It  follows  that  the  inverse  of  the 
control  effort  weighting  is  the  discrete  compensator  used  in 
the  feedback  loop.  While  any  discrete  compensator  can  be 
used,  the  PI  compensator  is  selected  because  its  charac¬ 
teristics  are  widely  known.  The  PI  compensator  used  is  of 
the  form: 

1  K  (  1  -z “  1  )  +  K. 

G  =  -  =  -2 - - - 1  (2.24) 

Q  (  1  -  z” 1 ) 

The  PI  compensator  is  used  instead  of  a  PID  compensator 
because  derivative  action  is  generally  not  required  for 
systems  without  time  delay.  As  the  self-tuning  controller 


' 
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predicts  the  output,  the  effect  of  process  time  delay  is 
removed  hence  derivative  action  is  not  necessary.  The 
incremental  PI  form  is  used  to  insure  that  Qu  tends  to  zero 

w 

at  steady  state  since  this  is  the  requirement  for  zero 
offset.  If  Qu  does  not  tend  to  zero  at  steady  state,  the 

V 

difference  between  the  weighted  setpoint  and  the  predicted 
output  ( Rw^_  -  (Pyt  +  k  )  )  is  different  from  zero  (cf.  Equation 
(2.23))  At  steady  state,  the  expected  predicted  output  is 
equal  to  the  actual  output.  If  Qu  does  not  equal  zero,  the 
difference  between  the  weighted  set  point  and  output  is 
different  from  zero  so  an  offset  occurs. 

The  control  law  can  now  be  formulated  in  terms  of  a 
predicted  scalar  output  defined  as: 

£5f+k  =  (Pyt  +  I<)*  "  Rwt  +  ^ut  (2.25) 

The  control  law  formulation  thus  involves  the  calculation  of 
u  by  zeroing  0^+k  at  every  sample  instant.  The  actual 
scalar  output  is  given  by: 


0 


t+k 


Rw  +  Qu 


(2.26) 


As  e.+k  and  Pyt+|<  are  uncor re  1  ated ,  the  error  in  the 
estimates  of  0^+k  can  be  obtained  from  Equation  (2.17)  as: 


0 


t+k 


-  £ 


t+k 


(2.27) 


Substitution  of  the  expression  for  ( Py  ^ )  from 
Equation  (2.16)  into  Equation  (2.25)  yields: 

C0*  =  F/P  y  +  EBu  +  CQu  -  CRw  +  ELzk‘kL  v  (2.28) 

t+k  d  t  t  t  t  i 


r! 


. 


. 
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Defining : 


y 


v 


k-kL 

Z  V 


H  =  -C 


D  =  EL 


(2.29) 

(2.30) 

(2.31 ) 

(2.32) 

(2.33) 


allows  the  control  law  to  be  written  as: 

C0^+^  =  Fy^.  +  EBu^  +  CQUj.  +  Hw^  +  Dv^  (2.34) 

Furthermore,  the  polynomial  associated  with  the  control 
ouput  u^_  can  be  defined  to  be  independent  of  the  control 
weight  function  by  setting: 


G  =  EB 


(2.35) 


which  then  allows  the  control  law  to  be  expressed  as: 

C0  £+  ^  =  Fy  ^  +  Gu^.  -  HQut  +  Hw^  +  Dv-t  (2.36) 

Alternatively,  all  the  polynomials  involving  u^_  can  be 
grouped  together  in  the  form  of: 


EB  +  QC 

G  =  - 

Q 


(2.37) 


So  the  control  law  takes  the  form: 

C0^+^  =  Fy^_  +  GQUj.  +  HWj.  +  Dvj. 


(2.38) 


. 
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The  block  diagram  of  the  control  law  as  expressed  in 
Equations  (2.36)  and  (2.38)  are  shown  in  Figures  2.2  and  2.3 
respect i ve 1 y . 


2.3.2  Closed  Loop  Response 

The  closed  loop  response  can  be  obtained  by 

* 

substituting  ( Pyt+k )  given  by  Equation  (2.17)  into  the 
control  law  Equation  (2.23).  After  rearranging  the  output, 
yt  is  given  by: 


y 


t 


1 

-  {EC  +  z_kRw.  -  z_i<Qu,} 

p  t  t  t 


(2.39) 


Solving  Equation  (2.1)  for  z'ku  and  substituting  into 

t 

Equation  (2.39)  yields  the  closed  loop  response: 


(EB  +  CQ)£t  +  z-kBRwt  +  z "  kLQLv  t 
(B  +  QA/P) 


(2.40) 


The  stability  of  the  system  is  governed  by  the 
char acter i s t i c  equation: 

PB  +  QA  =  0  (2.41 ) 

It  is  clear  that  if  B  is  nonminimum  phase;  i.e.  its  roots 
lie  outside  the  z-domain  unit  circle,  by  proper  choice  of  Q, 
the  roots  of  the  characteristic  equation  will  still  be 
inside  the  unit  circle.  The  closed  loop  response  is  then 
stable.  It  is  to  be  noted  that  as  in  the  case  of  the 
analytical  predictor  [63]  and  Smith  predictor  [64],  the 
effect  of  the  system  time  delay  (z”k)  has  been  removed  from 
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Figure  2.2  Block  Diagram  of  Self-Tuning  Controller  for  G=EB 
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Figure  2.3  Block  Diagram  of  Self-Tuning  Controller  for  G=(EB+CQ)/Q 
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the  characteristic  equation. 

2.3.3  Model  Following  For  Set  Point 

The  ability  of  the  algorithm  to  force  the  output  to 
follow  a  specified  response  can  be  clearly  seen  from  Figures 
2.2  and  2.3.  By  changing  the  weighting  functions  P  and  R, 
the  output  can  be  made  to  follow  a  desired  response. 
Furthermore,  the  set  point  weighting  function  R  can  be 
changed  without  altering  the  closed  loop  stability  of  the 
system.  In  contrast,  if  the  P  weight  is  changed,  it  will 
alter  the  character i st ic  equation  and  possibly  change  the 
stability  of  the  closed  loop  system. 

2.3.4  Controller  Parameter  Estimation 

The  algorithm  that  has  been  derived  is  for  a  system 
with  Known  structure  and  parameters.  Yet  in  practice,  this 
is  rarely  the  case.  Even  if  the  order  of  the  model  is 
Known,  usually  the  exact  process  dynamics  and  hence  the 
model  parameters  are  not  Known.  One  possible  solution  is  to 
use  an  on-line  identification  scheme  to  obtain  parameter 
estimates.  The  control  law  can  be  rendered  less  sensitive 
to  parameter  bias  and  variation  by  using  a  specific 
structured  '  control ler -model '  [24].  This  is  true  provided 

the  performance  function  is  minimized  when  the  estimated 
parameters  are  used.  Subtracting  Equation  (2.34)  from 
Equation  (2.27)  and  rearranging  yields  0t+k  in  linear 
regression  form,  that  is: 


I  V  IP 
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0 


t  +  k 


=  X 


th 


(1-00' 


t+k 


t  +  k 


where : 


(2.42) 


-t  =  f  Ft ’ Gt ■ Dt • Ht i 

h 

The  parameters  in  the  9  vector  can  now  be  estimated 
using  a  recursive  least  squares  algorithm  [65].  At  every 
sample  interval,  the  parameter  estimate  vector  0  is  updated 
according  to: 


0t  =  pyt  ■  Rwt-k  +  Qut-k 


e(  =  o  .  +  k .(0.  -  x‘  ,6.  o 

t  t-1  t  t  -t-k-t-1 


Ht 


=Ct-t-k 


!p  +  XTt.kPCtXt.k) 


D  c 

=  t+1 


-  £ct(' 


Hfct-k 


(2.43) 

(2.44) 


(2.45) 


(2.46) 


Here,  is  the  Kalman  gain  vector  and  P.c  is  a 
symmetric  matrix  that  is  proportional  to  the  covariance 
matrix.  To  allow  tracking  of  slow  time  varying  parameters, 
a  forgetting  factor  p  is  used  to  weight  out  the  effect  of 
past  data  or  alternatively,  a  positive  definite  matrix  can 
be  added  to  the  £c  update  relationship  to  weight  out  old 
data.  In  both  cases,  when  the  system  experiences  no 
disturbance  over  a  long  period,  the  £  matrix  can  become 
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extremely  large  [24].  This  occurs  because  without  any 
disturbances,  no  new  information  can  be  obtained  by  the 
identification  algorithm.  As  a  result,  the  Kalman  gains 
become  small  and  from  Equation  (2.46),  P  c  ( t )  approaches  the 
value  P_c(t-1)/p.  If  the  system  remains  undisturbed  over  a 
period  of  time  t,  then  £c ( t )  is  equal  to  Pc(to)/pt,  where  tQ 
refers  to  the  time  when  the  system  last  experienced  a 
disturbance.  As  the  forgetting  factor  (p)  is  less  than  one, 
as  time  passes,  1 / p ^  becomes  very  large  hence  causing  the 
matrix  to  'blow-up7  .  When  the  Pc  matrix  becomes  large, 
limited  machine  precision  may  cause  the  matrix  to  lose 
positive  definiteness.  As  a  result,  the  parameters  will 
diverge  causing  instability  in  the  control  loop.  The 
proposed  solution,  and  that  adopted  in  this  study,  is  to 
modify  the  identification  scheme  to  update  the  square  root 
of  the  Pc  matrix  [24].  The  square  root  will  always  be 
positive  definite  and  it  is  inherently  of  higher  numerical 
precision.  Consequently,  the  control  algorithm  now  consists 
of  using  a  recursive  square  root  identification  to  obtain 
the  parameter  vector  §  and  then  using  the  estimates  to 

V 

calculate  a  new  control  output  according  to  the  criteria 
( cf .  equation  (2.38)): 

xTtet  =  0  (2.471 

Astrom  has  shown  that  by  using  a  control  law  of  the 
form  given  by  Equation  (2.47),  all  vectors  with  the  form 
6  =  a0  wi 1 1  produced  the  same  control  [11].  This  means  that 
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the  parameters  will  not  necessarily  converge  to  provide  a 
unique  solution.  If  the  parameters  converge  to  an  a  that  is 
too  large  or  small,  numerical  problems  could  result.  In 
addition,  as  the  estimates  are  being  used  in  the  feedback 
controller,  they  are  correlated  with  the  disturbance.  It 
may  then  not  be  possible  to  identify  all  the  parameters  in 
the  0  vector.  Astrom  has  proposed  that  one  parameter  be 
fixed  to  overcome  these  problems.  One  obvious  choice  is  to 
fix  the  leading  coefficient  for  the  set  point  weighting 
polynomial  providing  w^  is  not  zero.  In  this  case,  hQ 
should  be  set  to  -1  as  it  is  equal  to  -c  .  Alternatively, 
Astrom  has  suggested  fixing  the  leading  coefficient  of  the  G 
polynomial;  i.e.,  the  scaling  parameter  for  the  controller 
output.  He  showed  that  if  this  parameter  is  much  smaller 
than  the  true  value,  the  other  parameters  will  diverge.  In 
contrast,  if  the  parameter  is  much  larger  than  the  true 
value,  the  control  action  becomes  very  small  and  convergence 
of  the  parameters  is  very  slow. 

2.3.5  Extension  to  Multivariable  case 

The  control  structure  derived  for  the  SISO  case  can  be 
extended  to  a  class  of  multivariable  systems  by  use  of  the 
feedforward  terms.  The  system  is  decoupled  by  using  feed¬ 
forward  action  of  the  manipulated  variable  of  one  output  as 
a  measurable  disturbance  for  the  other  outputs.  This 
treatment  will  be  restricted  to  systems  with  an  equal  number 
of  inputs  and  outputs.  In  addition,  the  time  delays  due  to 
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the  interactive  terms  must  be  greater  than  the  time  delay 
due  to  the  manipulated  variable  for  any  given  loop. 

Consider  a  multivariable  system  with  M-inputs  and 
M-outputs,  the  system  can  be  described  by  the  vector -matr i x 
difference  equation: 

A  8  .  L 

z'ki  J  t  A 1  Iyt  1  =  z'klj[|]  [uti  +  [C]  [£t]  +  z'kl'j[L][yt]  (2.48) 

where  [A],  [£] ,  [£] ,  [L]  are  polynomial  matrices  in  z"  ; 

[ytl  ,  [  u  t  ]  ,  [j^]  and  [yt]  are  vectors  correspondi  ng  to  the 

output,  input,  noise  and  feedforward  disturbances;  kA . . , 

^  J 

k8  •  •  .  k*-^  j  are  the  time  delays  of  the  output,  input  and 
feedforward  disturbance  for  the  jth  component  in  the  i th 
loop.  The  assumption  is  made  that  [B^]  is  nonsingular  and 
that  the  roots  of  all  the  row  polynomials  in  [C]  lie  inside 
the  z-domain  unit  circle. 

The  controller  performance  index  is  similar  to  the  SISO 

case : 


J  =  £{  UP][yt+k]  -  [£][wt])<[P][yt+k]  -  t  E 1  [  wt  ] )  T>  + 


E {  ((£'  I!ut])([fi'  ][ut])T} 


2.49) 


The  difference  is  that  ([P][yt  +  |J  "  [R][wt])  is  now  a  vector 
instead  of  a  scalar  value.  As  before,  [£] ,  [£]  ,  and  [Q'  ] 
are  matrix  weighting  functions  in  z”^.  Analogous  with  the 
SISO  case,  ( [£] [Yt+|J  )  can  be  replaced  by  the  predicted 
vector  (  [  P.  ]  [  ^  t+  k )  )  *  yielding  the  following  equation: 

J  =  E  {  ( ( [£][yt+kl)*  -  [£Hwt])  ( ( [£1  lyt+kl  *  *  ’  [EHwtl)Tt  + 

T  ,  I  It.  1 T  X  (2.50) 


E {  ([£'  ]  [u  ]  )  (  [£'  ][utl) 


+  let+k1  [^t+k!  T  } 


' 
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where  [e.  is  a  vector  of  weighted  prediction  errors  which 
is  assumed  to  be  uncorrelated  with  all  the  inputs  and 
outputs.  Minimizing  Equation  (2.50)  by  setting  3 J/3 I u t]  to 
zero  gives : 

{  U£]  [yt+k]  I  *  -  lii  [wt]}  [B  Q]  T  +  [S')[£'  ][ut]  =  0  (2.51) 

Defining: 

[£!  =  [£'  1  [Q'0]  (  [|0!  T)  -1  (2.52) 


the  control  law  can  now  be  written  in  terms  of  a  predicted 
vector  output  as: 


t  0*t+kJ  =  [£)  E  u  t  ]  +  {  U£Hyt+k])*  -  (£]  [  w  t  ] }  (2.53) 

Replacing  (  [  P.1  [  yt+ ^  ]  )  with  an  expression  analogous  to 
Eauation  (2.16)  and  using  similar  definition  to  those  in 
Equations  (2.29)  to  (2.33)  and  (2.37),  the  control  law 
becomes : 


- 1 


A 


B 


zM  j"ki  j  [G]  [ut  ]  =  - { z  k  i  j  -  k  i  j  [£] [y  ]  +  [H][wt]  +  [D][v  ]} 

m  ~  ifii  i/j  '  "  _t  =  _t 

(2.54) 

By  expressing  the  interactive  terms  as  feedforward 
disturbances,  the  parameters  can  be  identified  as  a  series 
of  M  loops.  Each  loop  would  correspond  to  the  parameters 
for  each  row  in  Equation  (2.54).  Using  recursive  least 
squares  identification,  the  algorithm  for  obtaining  the 
parameters  where  i  denotes  the  i th  loop,  is: 


P .  y . 


-  R  _.w 


i  it-kB 


+  Q  .u 


1  1 


i  t "  k  B  i  j 


(2.55) 


' 
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+  -u  (0it  '  Si't-i* 

(2.56) 

ii 

+-> 

^1 

=CitXit-kBii 

(2.57) 

(p,  +  £t1  b  £c,-  X-  B  1 

1  1 t- k? 1 t  1 t-kB. . 

11  11 

pc. 

=  H+i 

=  1/p.  *  (1  -  K.  X.  )£c • 

1  ’t  1 1 - k B  i  i  = 

(2.58) 

By  separating  the  loops,  the  identification  scheme  is 
simplified.  Futhermore,  the  weighting  functions  for  each 
loop  is  treated  separately,  that  is  the  control  effort 
weighting  Q.  can  be  tuned  individually  loop  by  loop.  In 
addition,  this  implementation  allows  the  used  of  a  different 
forgetting  factor  for  each  loop. 

For  most  processes,  the  control  law  given  by  Equation 
(2.54)  can  be  simplified  since  in  practice  one  output  in  a 
multivariable  system  does  not  directly  affect  another 
output.  It  is  the  controller  action  that  causes 
interactions.  It  follows  that  [F]  is  diagonal  and  the  terms 
kA  .  .  are  zero.  The  simplified  control  law  is  then: 

■  J 

[g][ut]  =  1/[fl]  {  [F][yt]  +  [H][wt]  +  [g][yt]}  (2.59) 

A  final  note,  although  the  identification  is  done  one 
loop  at  a  time,  in  general,  the  control  signals  have  to  be 
calculated  simultaneously.  The  [u  ]  vector  is  obtained 
directly  by  solving  Equation  (2.59).  However,  for  systems 
with  few  inputs  and  outputs,  the  solution  for  [Uj.1  can  be 
written  out  explicity  without  too  much  complication  so  that 
[u  ]  can  be  solved  by  direct  substitution. 
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3.  EXPERIMENTAL  IMPLEMENTATION 


3.1  DESCRIPTON  OF  EXPERIMENTAL  EQUIPMENT 

The  pilot  plant  distillation  column  in  the  Department 
of  Chemical  Engineering  at  the  University  of  Alberta  is  a 
22.5  cm  diameter  column.  It  has  a  thermosyphon  reboiler,  a 
total  condenser  and  eight  bubble  cap  trays  spaced  30.5  cm 
apart  with  four  bubble  caps  per  tray.  A  detailed 
description  of  the  column  and  the  associated  equipment  is 
given  by  Svrcek  [2]  and  Pacey  [6]. 

The  feed  is  a  water -methanol  mixture  containing  50%  by 
weight  methanol.  The  feed  enters  at  a  rate  of  18.08  g/s  and 
is  separated  to  an  overhead  composition  of  95%  by  weight 
methanol  and  a  bottom  composition  of  5%  by  weight  methanol. 
The  control  objective  is  to  minimize  variations  in  terminal 
compositions  when  the  feed  flow  rate  is  changed  by  25%  from 
its  normal  steady  state  feed  flow  rate. 

The  distillation  column  is  extensively  instrumented  to 
give  maximum  flexibility  for  studying  different  control 
stratergies.  A  schematic  diagram  showing  the  extent  of  the 
instrumentation  is  given  in  Appendix  A.  All  the  measurable 
variables  that  affect  the  performance  of  the  pilot  plant 
column  are  recorded.  Temperature  measurements  are  taken  for 
the  liquids  at  every  tray,  the  reboiler  and  condenser.  As 
well,  temperatures  and  flow  rates  of  the  bottoms, 
distillate,  feed,  reflux  and  steam  are  recorded.  The 
cooling  load  for  the  condenser  is  monitored  by  measuring  the 
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inlet  and  outlet  temperatures  and  flow  rate  of  the  cooling 
water.  The  overhead  composition  is  measured  continuously 
using  a  capacitance  probe  while  the  bottom  composition  is 
obtained  by  using  a  Hewlett-Packard  (HP)  model  5720A  gas 
chromatograph  equipped  with  an  automatic  sampling  valve. 

The  sampling  configuration  is  the  one  developed  by 
Bi lec  [9].  The  chromatogram  is  analyzed  on-line  by  a  HP 
1000  mi ni -computer  using  a  cycle  time  of  three  minutes. 
Additional  variables  that  are  measured  include  the  column 
pressure  and  the  liquid  levels  in  both  the  condenser  and 
reboiler.  All  the  signals  from  the  measured  variables  can 
be  brought  into  the  central  HP- 1000  mi ni -computer  located  in 
the  Data  Acquisition  Control  and  Simulinion  (DACS)  center  at 
the  Department  of  Chemical  Engineering  for  access  by  any 
program  executing  in  the  computer.  It  should  be  noted  that 
the  central  computer  is  not  the  same  one  used  for  analyzing 
the  gas  chromatograph.  Figure  3.1  shows  how  the 
distillation  column  is  interfaced  to  the  different  computers 
i n  the  DACS  center . 

The  column  is  instrumented  with  local  analog  control¬ 
lers  that  control  all  the  flow  rates,  the  condenser  and 
reboiler  levels,  the  feed  and  reflux  temperatures  and  the 
column  pressure.  The  instrumentation  is  such  that  control 
of  the  flows  and  temperatures  can  be  taken  over  by  the 
central  HP  1000  mi ni -computer  at  the  DACS  center.  The 
computer  can  assume  direct  digital  control  (DDC)  or  operate 
in  the  supervisory  mode.  With  DDC,  the  computer  sends  the 
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Figure  3.1  Schematic  Diagram  of  the  Distillation  Column  and  Interface 

to  Computers 
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control  output  to  the  valve  directly.  In  the  supervisory 
mode,  the  computer  relays  a  set  point  to  the  local  analog 
controller.  The  computer  provides  the  'master'  loop  and  the 
analog  controller  becomes  the  'slave'  loop. 

One  advantage  of  supervisory  control  over  DDC  is  that 
supervisory  control  can  remove  the  undesirable  effect  caused 
by  a  nonlinear  valve.  If  the  control  output  is  sent 
directly  to  a  nonlinear  valve,  then  the  same  change  in  valve 
position  for  different  inital  valve  position  will  result  in 
a  different  change  in  absolute  flow.  It  follows  that  for 
the  same  deviation  from  different  set  points,  the  change  in 
flow  caused  by  the  controller  will  be  different.  With 
supervisory  control,  the  appropriate  change  in  flow  is 
calculated  directly  in  the  computer  and  sent  to  the  flow 
controller.  The  required  change  in  valve  position  will  be 
subsequently  determined  by  the  analog  controller.  The 
result  is  that  the  same  output  will  always  produce  the  same 
change  in  the  flow. 

Another  advantage  with  supervisory  control  is  that  the 
flow  is  held  constant  in  between  sample  times  so  distur¬ 
bances  that  affect  the  flow  are  immediately  compensated  for 
by  the  analog  controller.  With  DDC,  the  control  output  is 
sent  directly  to  the  valve.  If  the  valve  position  is  held 
constant,  pressure  fluctuations  between  sample  times  can 
cause  the  flow  to  oscillate.  Local  control  is  advantageous 

in  the  event  of  computer  failure. 

However,  supervisory  control  becomes  less  desirable 
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when  the  response  time  of  the  slave  loop  is  of  the  same 
magnitude  as  the  sample  time  of  the  master  loop.  In  this 
situation,  the  dynamics  of  the  slave  loop  adds  directly  to 
that  of  the  master  loop  hence  causing  the  control  response 
to  deteriorate. 

For  the  pilot  plant  distillation  column,  the  response 
time  of  the  flow  loops  is  much  less  than  the  sample  time  of 
the  composition  loops.  Typically,  the  analog  controllers 
will  bring  the  flows  to  the  new  set  points  within  a  few 
seconds  which  is  significantly  less  than  the  sample  time  for 
the  composition  loops.  So  consequently,  the  supervisory 
mode  is  used  in  the  experimental  study. 

3.2  CONTROL  LAW  IMPLEMENTAION 

To  evaluate  the  self-tuning  controller  (STC),  the 
output  response  using  STC  will  be  compared  to  the  response 
achieved  using  a  well  tuned  PID  compensator.  The  comparison 
will  be  made  on  how  well  the  controllers  maintain  the 
terminal  compositions  when  the  feed  rate  is  changed. 

For  the  three  term  propor t i ona 1  - i ntegr a  1 -der i va t i ve 
compensator,  the  PID  option  in  the  Distributed  System  and 
Control  (DISCO)  package  is  used.  Details  of  the 
implementation  of  DISCO  on  the  pilot  plant  distillation 
column  and  the  PID  algorithm,  written  in  the  integral  form, 
are  given  by  Kan  [66].  The  STC  algorithm  could  not  be 
implemented  using  the  DISCO  executive  because  of  size 
restriction.  As  a  result,  a  completely  separate  program  had 
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to  be  written  to  execute  under  the  Real  Time  Executive 
operating  in  the  HP  1000  mi ni -computer .  A  listing  of  the 
program  is  given  in  Appendix  C.  The  STC  algorithm  used  in 
the  program  is  given  by  Equations  (2.55)  to  (2.59). 


4.  SELECTION  OF  MANIPULATED  VARIABLE 


4.1  INTRODUCTION 

Before  undertaking  the  control  study,  the  process 
output  and  control  input  pairs  for  the  pilot  plant 
distillation  column  were  first  chosen.  The  five  variables 
that  have  to  be  controlled  are:  top  composition,  bottom 
composition,  column  pressure,  level  in  the  reflux  drum  and 
level  in  the  reboiler.  The  five  available  manipulated 
variables  are:  reflux  flow,  distillate  flow,  bottom  product 
flow,  heat  input  to  the  reboi ler  and  heat  removal  by  the 
condenser.  The  variables  feed  temperature,  reflux 
temperature  and  feed  composition  will  be  kept  constant  for 
this  study. 

There  are  various  configurations  for  pairing  the 
controlled  and  manipulated  variables.  The  two  types  of 
strategies  that  have  been  used  most  frequently  are  the 
schemes  termed  material  balance  control  and  energy  balance 
control  [67].  With  the  energy  balance  scheme,  the  pressure 
of  the  column  is  controlled  by  manipulating  the  cooling  load 
in  the  condenser.  The  two  drum  levels  are  maintained  by 
changing  the  respective  output  flows.  The  top  composition 
is  controlled  by  the  reflux  flow  while  bottom  composition  is 
controlled  by  boilup  rate.  The  term  energy  balance 
originates  because  the  terminal  compositions  of  the  column 
are  controlled  by  changing  the  energy  balance  of  the  system 
and  it  is  clear  that  reflux  flow  and  reboiler  load  changes 
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will  directly  affect  the  energy  balance  of  the  column. 

With  material  balance  control,  the  terminal  composition 
is  controlled  by  directly  altering  the  material  balance  of 
the  distillation  column.  One  of  the  product  flow  rate  is 
manipulated  to  control  its  own  composition.  Level  for  that 
end  of  the  column  is  then  controlled  by  either  reflux  flow 
or  reboiler  load  depending  on  which  product  flow  rate  is 
used.  The  remaining  variables  will  be  paired  the  same  way 
they  are  matched  in  the  heat  balance  stratery. 

Distillate  flow  is  usually  selected  to  control  the 
material  balance.  The  reason  for  this  is  that  problems 
could  develop  if  the  the  reboiler  level  is  to  be  controlled 
by  manipulation  of  the  reboiler  load.  As  the  distillate 
flow  is  usually  used,  the  ensuing  discussion  on  material 
balance  control  will  be  based  on  using  distillate  flow  to 
control  the  top  composition.  Analogous  arguments  can  be 
given  that  will  be  applicable  to  the  case  when  the  bottom 
composition  is  controlled  by  manipulating  bottom  product 
flow  rate.  Schematic  diagrams  of  the  overhead  composition 
control  configuration  for  the  material  and  energy  balance 
control  schemes  are  shown  in  Figures  4.1  and  4.2. 

The  concept  of  'floating'  the  reflux  flow  on  level 
control  is  not  an  easy  one  to  accept  since  basic  to  the 
distillation  column  design  procedure  is  the  assumption  of 
the  operation  at  a  constant  reflux  flow  and  fixed  reflux 
ratio.  Nevertheless,  studies  on  different  distillation 
columns  have  shown  that  the  use  of  the  material  balance 
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Figure  4.1  Schematic  Diagram  of  the  Energy  Balance  Control  Scheme 
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Figure  4.2  Schematic  Diagram  of  the  Material  Balance  Control  Scheme 
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control  scheme  can  result  in  improved  control  behaviour  over 
the  case  when  energy  balance  control  is  used  [68-73].  In 
certain  cases,  it  has  been  found  that  it  is  not  possible  to 
stabilize  control  without  using  the  material  balance 
approach  [72,73].  In  both  cases,  control  of  the  overhead 
composition  could  not  be  realized  by  direct  manipulation  of 
reflux  flow  because  of  the  very  high  external  reflux  ratio 
( greater  than  25 ) . 

The  main  advantage  of  the  material  balance  strategy  is 
that  it  exhibits  inherent  reflux  control  when  an  energy 
disturbance  is  encountered  because  the  reflux  flow  readjusts 
itself  before  any  action  is  taken  by  the  composition 
controller.  This  is  accomplished  because  an  energy 
disturbance  will  change  the  vapor  rate  in  the  column.  As 
the  new  vapor  rate  reaches  the  condenser,  the  reflux  drum 
level  starts  to  change.  Since  the  reflux  flow  is  used  to 
control  the  level,  it  will  be  adjusted  to  maintain  level  set 
point.  The  overall  result  is  that  reflux  flow  will  change 
in  the  same  direction  as  the  change  in  vapor  rate.  Since 
the  scheme  inherently  handles  heat  disturbances,  material 
balance  control  reduces  the  interaction  that  occurs  during 
simultaneous  control  of  top  and  bottom  composition.  Varia¬ 
tions  in  the  reboiler  load  due  to  the  bottom  composition 
controller  are  treated  as  heat  disturbances  by  the  overhead 
controller  so  the  reflux  flow  will  be  adjusted  in  the 
appropriate  direction  before  any  change  in  overhead 
composition  is  detected. 
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Other  advantages  of  the  material  balance  scheme  are  the 
fol lowing: 

a.  For  the  case  where  the  reflux  ratio  is  very  large, 
using  the  distillate  flow  to  control  the  reflux 
drum  level  can  be  very  difficult.  Since  the 
distillate  flow  rate  is  much  lower  than  the  reflux 
flow,  a  large  percentage  change  in  distillate  flow 
is  required  to  compensate  for  small  percentage 
changes  in  reflux  flow.  These  large  swings  in 
distillate  flow  can  be  a  big  problem  especially  in 
situations  where  the  distillate  product  is  a  feed 
to  another  column.  Additional  drum  capacity  might 
be  needed  to  damp  out  the  variation. 

b.  The  material  balance  stratery  is  readily  adapted 
for  feedforward  application  since  the  required 
distillate  to  feed  ratio  can  be  determined  exactly 
from  the  material  balance  relationship.  Thus,  the 
necessary  change  in  distillate  flow  to  compensate 
for  changes  in  feed  flow  rate  and  feed  composition 
can  be  easily  calculated  and  subsequently  fed 
forward  through  the  appropriate  dynamic 
compensator . 

c.  The  material  balance  scheme  will  not  exhibit  a 
transient  inverse  response  that  is  sometimes 
observed  when  reflux  flow  is  directly  used  to 
control  overhead  composition.  The  inverse  reponse 
as  explained  by  RosenbrocK  [74]  is  caused  by  plate 
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overflow  when  an  increased  vapor  rate  first 
reaches  the  top  of  the  distillation  column.  The 
overflow  acts  as  an  increased  reflux  flow  hence 
causing  the  composition  to  initially  become  richer 
in  the  more  volatile  component.  When  the  overflow 
stops,  the  increased  vapor  rate  takes  over  and  the 
composition  becomes  leaner  in  the  more  volatile 
component.  With  material  balance  control,  the 
increased  vapor  rate  would  cause  the  reflux  to 
increase  and  as  a  result,  the  inverse  response 
should  not  occur. 

The  main  disadvantage  of  the  material  balance  strategy 
over  the  heat  balance  control  case  is  the  addition  of  the 
control  dynamics  of  the  level  loop  into  the  dynamics  of  the 
composition  loop.  When  the  top  composition  varies  from  the 
set  point,  the  distillate  flow  is  first  adjusted.  This 
change  in  turn  affects  the  level  in  the  reflux  drum. 
Subsequently,  the  reflux  flow  then  adjusts  to  maintain  the 
drum  level.  The  composition  will  not  begin  to  change  until 
the  reflux  flow  has  been  changed.  The  overall  result  is 
that  more  dynamics  is  added  hence  slowing  the  response  of 
the  overhead  composition  loop.  Futhermore,  controller 
response  to  set  point  changes  could  be  slower  than  the  case 
when  the  reflux  is  adjusted  directly. 


46 


4.2  ANALYSIS  OF  INTERNAL  REFLUX  RATIO  SENSITIVITY  TO  HEAT 
DISTURBANCE 

The  sel f -regul at  ion  of  the  reflux  flow  in  material 
balance  control  will  lead  to  improved  control  performance 
only  if  the  change  in  reflux  flow  is  changed  in  a  way  which 
maintains  the  internal  reflux  ratio.  Clearly,  it  is  the 
maintenance  of  the  internal  reflux  ratio  that  allow  the 
material  balance  scheme  to  be  less  sensitive  to  heat 
disturbances.  By  holding  a  constant  internal  reflux  ratio, 
the  slope  of  the  operating  curve  remains  unchanged  hence  the 
composition  is  unaltered. 

For  a  distillation  column  operating  with  a  large 
external  reflux  ratio,  a  material  balance  strategy  will 
maintain  the  internal  reflux  ratio  better  than  the  energy 
balance  case.  However,  if  the  external  reflux  ratio  is  low, 
the  sel f-regulat ing  change  in  reflux  flow  can  cause  the 
internal  reflux  ratio  to  vary  more  than  the  energy  balance 
scheme.  The  best  way  to  illustrate  the  effect  is  a  simple 
numerical  example.  It  should  be  noted  that  in  the  ensuing 
analysis,  the  V/L  ratio  will  be  considered  instead  of  L/V 
traditionally  defined  as  the  internal  reflux  ratio.  The 
reason  for  this  is  that  the  region  of  interest  is  for 
internal  reflux  ratio  less  than  1.  The  reciprocal,  V/L, 
will  have  a  greater  numerical  range. 

Consider  a  situation  where  the  external  reflux  ratio  is 
10.0,  that  is: 

Reflux  Flow  (L) 


Distillate  Flow  ( D ) 


10.0 


(4.  1  ) 
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or  , 


L  =  10. CD 


(4.2) 


Since  the  vapor  rate  (V)  is  given  by: 


V  =  L  +  D 


(4.3) 


Combining  Equation  (4.2)  and  (4.3)  yields: 


V  =  1 1 . OD 


(4.4) 


and  the  resulting  ratio  is: 


V  11. OD 


1  .  1 


(4.5) 


L  10. OD 

Suppose  that  the  column  is  subjected  to  energy 
disturbances  that  correspond  to  a  required  change  in  the 
vapor  rate  of  a  20%  increase  and  a  20%  decrease.  The 
required  new  V/L  ratio  for  both  material  balance  and  energy 
balance  stratergies  for  such  changes  are  given  in  Table  4.1. 
Recall  that  for  material  balance  control,  the  reflux  flow 
will  change  by  the  same  amount  the  vapor  rate  is  changed. 
From  Equation  (4.4),  a  20%  change  in  vapor  rate  is: 


(4.6) 


0.2V  =  0.2  *  1 1 .0D  =  2. 2D 


It  is  evident  from  the  values  in  Table  4.1  that  the 
material  balance  control  clearly  maintains  a  more  constant 
V/L  ratio  than  the  heat  balance  case  as  the  ratio  under  the 
material  scheme  changes  less  than  3%  compared  to  the  20% 


48 


changes  for  the  energy  balance  case. 


Table  4.1  Variation  of  V/L  with  Control  Scheme  for  ±20% 
Changes  in  Vapor  Rate  at  an  External  Reflux 
Ratio  of  10.0 


Control 

Pos i t i ve 

Deviation 

Scheme 

Disturbance 

From  S.S. 

Mater i a  1 

11.0+2.2 

-  =1.08 

-1.6% 

Ba 1 ance 

10.0+2.2 

Energy 

11.0+2.2 

=  1.32 

+20.0% 

Ba 1 ance 

10.0 

Negative 
Di sturbance 


1 1 .0-2.2 

-  =  1.13 

10.0-2.2 


11.0-2.2 

-  =0.88 

10.0 


Devi  at i on 
From  S.S 


+  2.6% 


-20.0% 


Now,  consider  the  case  when  the  external  reflux  ratio 
is  only  1.0,  that  i s : 

-  =  1.0  (4.7) 

D  • 


so , 


V  =  L  +  D  =  2.0D 


(4.8) 


and  the  steady  value  of  V/L  is: 


V  2.0D 


L 


D 


2.0 


(4.9) 


■ 


. 
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The  new  V/L  ratios  that  result  from  ±20%  changes  in 
vapor  rate  are  given  in  Table  4.2. 


Table  4.2  Variation  of  V/L  with  Control  Scheme  for  ±20% 
Changes  in  Vapor  Rate  at  an  External  Reflux 
Ratio  of  1.0 


Control 

Scheme 


Posi t i ve 
Di sturbance 


Deviation  Negative  Deviation 

From  S.S.  Disturbance  From  S.S 


Material  2. 0+0.4 

-  =  1.71 

Ba 1 ance  1  .  0  +  0 . 4 


2. 0-0.4 

-14.3%  -  =  2.67  +33.3% 

1 .0-0.4 


Energy  2. 0+0.4 

-  =2.40 

Ba 1 ance  1  .  0 


2. 0-0. 4 

+  20.0%  -  =  1.60  -20.0% 

1.0 


As  can  be  seen  from  the  numerical  values  in  Table  4.2, 
which  are  markedly  different  from  those  in  Table  4.1,  the 
material  balance  control  scheme  is  preferable  for  a  postive 
disturbance  but  just  the  reverse  for  the  negative 
disturbance.  Cleary,  for  a  situation  where  the  external 
reflux  ratio  is  1.0,  the  energy  balance  strategy  will  prove 
to  be  less  sensitive  to  heat  disturbances.  The  choice  of 
energy  balance  control  becomes  even  more  advantageous  if  it 
is  recalled  that  with  the  energy  balance  scheme,  the 
dynamics  of  the  reflux  drum  level  do  not  enter  the  overhead 
composition  control  loop. 

These  simple  numerical  calculations  have  clearly 
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illustrated  that  the  material  balance  control  scheme  will 
not  always  be  superior  to  heat  balance  control  with  the 
choice  of  scheme  dependent  on  the  external  reflux  ratio. 
While  it  is  true  that  most  industrial  columns  operate  with 
an  external  ratio  much  greater  than  one,  that  is  not 
necessarily  the  case  with  pilot  plant  distillation  columns. 
In  this  study,  the  pilot  plant  column  operating  conditions 
are : 


Reflux  Flow  =  10.2  g/s 
Distillate  Flow  =  8.8  g/s 


so,  the  steady  state  V/L  ratio  is: 


V  10.2  +  8.8 

-  =  -  =1.87 

L  10.2 


The  change  in  the  V/L  ratio  for  ±20%  changes  in  the 
vapor  rate  are  shown  in  Table  4.3. 


Table  4.3  Variation  of  V/L  with  Control  Scheme  for  ±20% 
Changes  in  Vapor  Rate 


Control 

Scheme 


Posi t i ve 
Di sturbance 


Deviation  Negative  Deviation 
From  S.S.  Disturbance  From  S.S. 


Material  19.0+3.8 

- =1.63 

Balance  10.2+3.8 


19.0+3.8 

-12.6%  - =2.34  +27.6% 

10.2-3.8 


Energy  19.0+3.8 

- =2.24 

Balance  10.2 


19.0-3.8 

+  20.0%  - =1.49  -20.0% 

10.2 
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The  results  presented  in  Table  4.3  indicate  that  the 
reduced  deviation  with  the  material  balance  control  compared 
to  energy  balance  control  for  a  positive  disturbance  (7.4% 
lower  deviation)  is  just  about  equal  to  the  higher  deviation 
of  7.6%  for  the  negative  disturbance.  As  both  stratergies 
are  equally  sensitive  to  heat  disturbances,  heat  balance 
control  is  the  strategy  chosen  because  the  material  balance 
scheme  has  the  added  disadvantage  of  having  a  slower 
response  due  to  the  inclusion  of  the  condenser  level 
dynamics . 


4.3  EXPERIMENTAL  OPEN  LOOP  RESPONSES 

To  substantiate  the  analysis  carried  out  in  the 
previous  section,  experimental  open  loop  responses  where 
obtained  for  the  top  and  bottom  compositions  to  changes  in 
steam  flow  using  both  the  material  balance  and  energy 
balance  configurations.  The  steam  flow  changes  simulated 
heat  disturbances  entering  the  system. 

For  each  control  strategy,  steam  flow  rate  was 
increased  and  decreased  to  flow  rates  representing  a  change 
of  7  1/2%  from  its  steady  state  value.  The  open  loop 
responses  for  the  material  balance  scheme  are  shown  in 
Figures  4.3  and  4.4.  The  corresponding  responses  for  energy 
balance  control  are  shown  in  Figures  4.5  and  4.6.  The 
results  are  summarized  in  Table  4.4. 
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Table  4.4  Top  and  Bottom  Composition  Changes  to  Steam  Flow 
Rate  Disturbances. 


7  1/2%  Increase  In  7  1/2%  Decrease  in 

Steam  Flow  Rate  Steam  Flow  Rate 


Control 

XD 

XB 

XD 

XB 

Scheme 

Change 

Change 

Change 

Change 

Mater i a  1 

Ba 1 ance 

+0.82% 

-0.55% 

-  1  . 32% 

+1.32% 

Energy 

Ba 1 ance 

-1  .50% 

-4.08% 

+0.58% 

+6.57% 

The  results  of  Table  4.4  are  in  agreement  with  the 
predictions  given  by  the  analysis  of  the  vapor- to- 1 iquid 
ratio  for  vapor  rate  changes.  The  values  in  Table  4.3 
indicate  that  the  overhead  composition  should  be  more 
sensitive  to  a  positive  energy  disturbance  when  the  energy 
balance  scheme  is  used  than  if  material  balance  control  is 
used.  This  is  indeed  the  case  as  can  be  seen  from  the  1.5% 
decrease  in  overhead  composition  with  energy  balance  control 
compared  to  the  0.82%  increase  under  the  material  balance 
scheme.  The  results  for  the  negative  disturbance  are  also 
consistent  with  those  in  Table  4.3  since  the  top  composition 
decreases  by  1.32%  for  material  balance  control  but  only 
increases  by  0.58%  for  energy  balance  control. 

The  potential  advantage  of  the  material  balance  scheme 
for  systems  with  a  large  external  reflux  ratio  is  illus¬ 
trated  by  considering  the  bottom  composition.  Since  the 
entering  feed  is  a  liquid,  the  liquid  rate  below  the  feed 
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Figure  4.3  Open  Loop  Response  to  a  7  1/2%  Step  Decrease  in 

Steam  Flow  With  A  Material  Balance  Configuration 
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Figure  4.4  Open  Loop  Response  to  a  7  1/2%  Step  Increase  in 

Steam  Flow  With  A  Material  Balance  Configuration 
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Figure  4.6  Open  Loop  Response  to  a  7  1/2%  Step  Increase  in 

Steam  Flow  With  an  Energy  Balance  Configuration 
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plate  is  increased  substantially.  In  effect,  more  reflux  is 
being  added  to  the  column.  From  the  results  of  Table  4.1, 
it  should  follow  that  the  vapor  -  to- 1 iquid  ratio  and  hence 
the  bottom  composition  should  be  less  sensitive  to  heat 
disturbances  if  material  balance  control  is  used.  This  is 
in  fact  the  situation  observed.  For  both  positive  and 
negative  disturbances,  there  is  a  more  significant  change  in 
the  bottom  composition  with  the  energy  balance  configuration 
than  with  material  balance  control.  For  the  positive 
disturbance,  the  composition  change  with  energy  balance 
control  is  7  1/2  times  larger  than  the  change  with  material 
ba 1 ance  control . 


4.4  CONCLUSION 

Material  balance  control  is  shown  to  be  superior  to 
energy  balance  control  for  the  case  of  a  large  external 
reflux  ratio  while  with  smaller  reflux  ratio  it  will  be  more 
advantageous  to  use  the  energy  balance  scheme.  Comparison 
of  the  change  in  the  vapor- to- 1 iquid  ratio  for  the  same 
change  in  vapor  rate  for  both  control  schemes  can  help  to 
decide  which  strategy  should  be  used  for  a  given  situation. 

For  the  column  operating  conditions,  analysis  of  the 
vapor- to- 1 iquid  ratio  for  vapor  rate  changes  and  subsequent 
experimental  verification  have  shown  that  the  choice  of 
control  strategy  depends  on  the  type  of  disturbance.  Energy 
balance  control  is  less  sensitive  to  a  negative  energy 
disturbance  but  because  the  material  balance  control  scheme 
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inherently  includes  condenser  level  dynamics,  the  energy 
balance  scheme  is  employed  in  this  study. 


5.  EXPERIMENTAL  RESULTS 


5.1  INTRODUCTION 

This  chapter  presents  the  results  of  the  experimental 
runs  that  have  been  carried  out  on  the  distillation  column 
using  the  PID  and  self-tuning  controllers.  A  brief 
discussion  of  the  experimental  procedures  and  the  actual 
programs  that  were  used  will  be  given.  The  results  are 
presented  in  the  following  order: 

a.  PID  runs  for  top  composition,  bottom  composition 
and  dual  composition  control; 

b.  STC  runs  for  top  composition  control; 

c.  STC  runs  for  bottom  composition  control; 

d.  STC  runs  for  dual  composition  using  multi  loop 
configuration  (ML-STC); 

e.  STC  runs  for  dual  composition  using  multivariable 
configuration  (MV-STC). 

Control  of  top  product  composition  will  be  used  to 
study  most  of  the  properties  of  the  self-tuning  controller 
while  bottom  composition  will  be  used  to  evaluate  the 
feedforward  (FF)  feature  of  the  STC. 

5.2  EXPERIMENTAL  PROCEDURE 

Four  types  of  disturbances  were  used  to  perturb  the 

distillation  column: 

a.  a  25%  reduction  in  feed  flow  rate; 

b.  a  subsequent  feed  flow  rate  increase  to  return  the 
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rate  to  its  normal  steady  state  flow; 

c.  a  25%  increase  in  feed  flow  rate; 

d.  a  subsequent  feed  flow  rate  decrease  to  return  the 
rate  to  its  normal  steady  state  value. 

For  the  study  of  top  composition  control,  the 
disturbances  of  a  25%  increase  in  feed  flow  rate  and  the 
subsequent  return  of  the  flow  rate  to  its  steady  state  value 
were  not  used  because  the  top  composition  showed  virtualy  no 
change  for  these  disturbances  as  can  be  seen  from  the 
responses  in  Figures  5.1  and  5.2  (The  arrow  in  the  figures 
indicate  the  time  that  the  step  change  in  feed  flow  rate  is 
introduced).  In  addition,  only  a  limited  number  of 
multi  loop  and  multivariable  runs  were  performed  for  these 
two  disturbances  because  a  condition  of  flooding  nearly 
resulted  when  the  feed  rate  was  increased  25%  above  the 
steady  state  flow  rate. 

With  single  loop  control,  the  self-tuning  controller 
was  started  with  zero  initial  parameter  estimates  and  the 
initial  £c  matrix  is  set  at  100CH  to  show  that  there  was 
lack  of  confidence  in  the  initial  estimates.  The  controller 
was  given  time  to  obtain  reasonable  estimates  before  a  feed 
disturbance  was  introduced.  For  the  multi  loop  and 
mu ti variable  configuration,  the  self-tuning  controller  was 
started  with  the  final  estimates  obtained  from  the  single 
loop  control  tests.  To  show  that  there  was  confidence  in 
these  estimates,  the  initial  P_c  matrix  was  set  equal  to 
.OH.  As  with  the  single  loop  control  case,  the  parameters 
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Figure  5.1  PID  Control  of  Top  Composition  for  a  +25%  Step 

Change  in  Feed  Flow  Rate 
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were  allowed  to  converge  before  a  feed  change  was  made. 

When  comparing  the  control  performance  of  the  self- tuning 
controller  with  that  of  the  PID  controller,  only  STC  runs 
for  which  the  disturbance  is  introduced  after  the  initial 
identification  are  considered 

The  discussion  given  in  Section  2.3.1  has  shown  that 
the  self-tuning  controller  can  be  thought  of  as  a  feedback 
controller  with  a  discrete  compensator  which  uses  the 
predicted  output  instead  of  the  measured  output  as  the 
feedback  signal  for  control  action.  From  Equation  (2.23), 
it  can  be  seen  that  the  discrete  compensator  is  equal  to  the 
inverse  of  the  Q-weighting  (Q-WT).  By  selecting  the  inverse 
of  the  Q-weighting  to  be  the  PI  controller  and  assuming  the 
predicted  output  to  be  correct,  standard  techniques  proposed 
by  Verbruggen  et  al.  [75]  can  be  used  to  obtain  an  initial 
choice  for  the  PI  constants  based  on  the  open  loop  response. 
Since  the  predicted  output  is  used,  the  PI  constants  are 
those  based  on  the  system  except  for  the  effect  of  the 
process  time  delay. 

For  the  single  loop  control  tests,  the  constants 
obtained  from  the  method  based  on  the  open  loop  responses 
produced  good  control  but  it  was  found  that  increasing  the 
integral  gain  improved  the  control  response.  For  the 
multi  loop  and  multivariable  cases,  a  certain  amount  of 
de-tuning  was  necessary  to  prevent  oscillatory  responses. 

Experience  with  tuning  the  Q-weighting  showed  that  the 
response  was  not  as  sensitive  to  variation  in  the 
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Q-weighting  parameters  compared  to  the  sensitivity  of  the 
PID  compensator  to  changes  in  the  PID  constants.  The  reason 
for  this  insensitivity  is  due  to  the  fact  that  the 
parameters  change  to  compensate  to  some  degree  for  changes 
in  the  Q-weighting  constants.  Although  the  sensitivity  was 
different,  it  was  found  that  the  PI  constants  in  the 
Q-weighting  could  be  tuned  in  the  same  way  a  normal  PI 
compensator  is  tuned;  i .e. ,  the  relationship  between  the 
proportional  and  integral  gains  and  how  they  affect  the 
output  response  remains  the  same. 


5.3  COMPUTER  PROGRAMS 

With  the  change  of  computer  system  in  the  DACS  center, 
support  programs  for  the  distillation  column  had  to  be 
written  for  the  new  system.  A  monitor  program,  BDC99,  was 
written  to  monitor  Key  process  variables  on  the  column. 

When  the  values  of  the  variables  exceeds  specified 
constraints,  BDC99  will  shut  down  the  distillation  column. 

As  discussed  in  chapter  three,  the  chromatogram 
analysis  is  done  by  a  computer  seperate  from  the  central 
computer  where  the  other  process  measurements  are  available. 
As  the  control  algorithm  will  execute  in  the  central 
computer,  a  program,  GCLNK  had  to  be  written  to  receive  the 
report  from  the  GC  computer  and  transfer  the  information 
into  some  buffer  in  the  central  computer.  Programs  GCSCH 
and  GCSWT  then  processed  the  information  in  the  buffer  to 
extract  the  value  of  the  the  bottom  composition  from  the 
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report.  The  programs  also  served  to  start  the  GC  sample 
cycle  at  the  appropriate  time. 

The  implementation  of  the  self-tuning  control  algorithm 
on  the  column  was  carried  out  with  programs  STC  and  HSET. 
Program  HSET  is  used  to  alter  the  data  file  that  contains 
the  starting  parameters  and  options  used  for  the  self-tuning 
controller  while  program  STC  contains  the  control  algorithm. 
In  addition  to  the  algorithm  itself,  program  STC  takes  care 
of  storing  the  relevant  variables  in  a  disk  file.  When  the 
STC  program  is  executing,  communication  with  the  program  can 
be  accomplished  through  a  group  of  digital  switches. 

Utility  programs  HPLOT  and  ERSUM  were  written  to  provide 
plots  and  calculation  of  the  sum  of  the  absolute  error  (SAE) 
of  the  results  stored  in  the  data  files  for  a  user  specified 
period  of  time  (number  of  sampling  intervals). 

A  listing  of  programs  HSET  and  STC  is  included  in 
Appendix  C.  A  listing  of  the  remaining  programs  is  only 
included  in  certain  copies  of  this  work. 


5.4  CONTROL  BEHAVIOUR  USING  PID  CONTROL 

The  control  performance  obtained  using  conventional  PID 
control  action  will  be  used  as  the  base  case  for  evaluation 
of  the  self-tuning  controller.  For  the  three  cases  of:  top 
composition  control,  bottom  composition  control  and 
simultaneous  product  composition  control,  the  PID  results 
presented  are  for  a  well  tuned  PID  compensator.  While  the 
performance  cannot  be  guaranteed  to  have  been  obtained  using 


66 


true  optimal  tuned  PID  controller  settings,  the  constants 
have  been  so  well  tuned  that  any  changes  in  their  values 
cause  the  control  performance  to  deteriorate.  The 
controller  performance  is  measured  by  the  sum  of  the 
absolute  error  (SAE).  For  top  composition  control,  with  a 
sampling  time  of  1  minute,  the  SAE  is  taken  over  one  hour  or 
60  samples  after  the  feed  flow  rate  change.  For  bottom  and 
dual  composition  control,  with  a  sampling  time  of  3  minute, 
the  SAE  is  taken  over  1  1/2  hour  period  or  30  samples  after 
the  feed  flow  rate  change. 

For  the  single  loop  control  test,  the  first  choice  of 
PID  constants  was  obtained  using  the  techniques  proposed  by 
Verbruggen  et.  al.  [75]  which  bases  the  selection  of  PID 
constants  on  the  system  open  loop  response.  Using  these 
values  as  a  starting  point,  the  parameters  were  altered  in  a 
systematic  way  until  further  changes  in  the  constants 
resulted  in  deter iorat ion  of  the  control  performance.  For 
the  simultaneous  control  case,  the  well  tuned  single  loop 
constants,  reduced  by  one  half  were  used  as  the  initial 
choice.  The  constants  were  reduced  because  of  the 
interaction  between  the  top  and  bottom  composition  control 
loops . 

The  experimental  results  for  top  composition  control 
are  those  given  in  Section  5.2  (cf.  Figures  5.1  and  5.2)  and 
in  Figures  5.3  and  5.4.  The  bottom  composition  control 
results  are  shown  in  Figures  5.5  to  5.8.  For  bottom 
composition  control,  the  PID  constants  are  significantly 
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Figure  5.3  PID  Control  of  Top  Composition  for  a  -25%  Step 

Change  in  Feed  Flow  Rate 
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Figure  5.4  PID  Control  of  Top  Composition  for  a  Step 

Increase  in  Feed  Flow  Rate  to  its  Normal  Steady 
State  Value 
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PID  Control  of  Bottom  Composition  for  a  +25% 
Step  Change  in  Feed  Flow  Rate 
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Figure  5.6  PID  Control  of  Bottom  Composition  for  a  Step 

Decrease  in  Feed  Flow  Rate  to  its  Normal  Steady 
State  Value 
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Figure  5.7  PID  Control  of  Bottom  Composition  for  a  -25% 

Step  Change  in  Feed  Flow  Rate 
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Figure  5.8  PID  Control  of  Bottom  Composition  for  a  Step 

Increase  in  Feed  Flow  Rate  to  its  Normal  Steady 
State  Value 


} 


73 


different  for  the  positive  and  the  negative  feed  flow  rate 
disturbances  as  can  be  seen  from  Table  5. 1 .  This  is  due  to 
the  nonlinear  dynamic  characteristics  of  the  column.  The 
results  for  the  dual  composition  control  are  shown  in 
Figures  5.9  to  5.12.  The  PID  constants  are  tuned  on  the 
basis  of  a  -25%  step  change  in  feed  flow  rate  only.  A 
summary  of  the  SAE  for  these  results  is  given  in  Table  5.1. 


Table  5.1  Summary  of  Absolute  Error  Values  For  PID  Control 


Type  of 
Control 

Type  of 

D i sturbance 

F igure 

SAE 

(wt%) 

PID  Constants 
kp  /ki  /kd 

Top 

+25%  to  feed 

5 .  1 

0.460 

2.81/0.73/0.65 

Top 

Return  to  S.S. 

5.2 

0.785 

2.81/0.73/0.65 

Top 

-25%  to  feed 

5.3 

3.870 

2.81/0.73/0.65 

Top 

Return  to  S.S. 

5.4 

3.410 

2.81/0.73/0.65 

Bot  tom 

+25%  to  feed 

5.5 

20.72 

-.50/-. 13/-. 08 

Bot  tom 

Return  to  S.S. 

5.6 

26.05 

-.50/-. 13/-. 08 

Bot  tom 

-25%  to  feed 

5.7 

18.43 

-.70/-. 15/-. 67 

Bot  tom 

Return  to  S.S. 

5.8 

14.00 

-.70/-. 15/-. 67 

Dua  1 

+25%  to  feed 

5.9 

top  /  bot 
2.75/27.3 

top  constants: 

Dua  1 

Return  to  S.S. 

5.  10 

2.87/27.1 

3.57/0.93/0.83 

Dua  1 

-25%  to  feed 

5.11 

2.13/21.7 

bot  constants: 

Dua  1 

Return  to  S.S. 

5.  12 

2.30/23.0 

-.58/-.  13/-.  13 
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Figure  5.9  Multi  loop  PID  Control  of  Terminal  Compositions 

for  a  +25%  Step  Change  in  Feed  Flow  Rate 
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Figure  5.10  Multi  loop  PID  Control  of  Terminal  Compositions 

for  a  Step  Decrease  in  Feed  Flow  Rate  to  its 
Normal  Steady  State  Value 
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Figure  5.11  Multi  loop  PID  Control  of  Terminal  Compositions 

for  a  -25%  Step  Change  in  Feed  Flow  Rate 
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5.5  TOP  COMPOSITION  CONTROL 

This  section  contains  the  experimental  responses 
achieved  using  the  self-tuning  controller  for  overhead 
composition  control.  Results  are  presented  to  show  the 
effect  on  control  behaviour  of: 

a.  removing  Q-weighting  (Q-WT); 

b.  changing  sampling  time; 

c.  stopping  identification. 

Top  composition  control  was  selected  to  demonstrate 
some  characteristic  of  the  self-tuning  control ler  because  of 
the  loop's  smaller  time  constant  and  the  allowable  user 
choice  of  sampling  rate  for  the  overhead  composition  from 
the  continuos  measurement  signal  of  the  capacitance  cell. 

5.5.1  Control  Performance  Using  the  STC  with  Q-weighting 

The  initial  controller  parameters  and  options  used  for 
self- tuning  control  when  the  column  was  subjected  to  the 
different  feed  disturbances  are  shown  in  Table  5.2. 

The  top  composition  response  when  the  self- tuning 
controller  is  first  made  operative  with  zero  initial 
parameter  estimates  is  shown  in  Figure  5.13.  The  changes  in 
the  parameters  during  the  initial  tuning  period  are  shown  in 
Figures  5.14  and  5.15.  It  can  be  seen  that  even  starting 
with  zero  inital  estimates,  the  parameters  have  converged 
after  less  than  30  sampling  intervals.  With  the  parameters 
having  converged,  the  feed  flow  rate  was  decreased  by  25%, 
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Table  5.2  Controller  Parameters  and  Options  for  STC  Top 
Composition  Control 


Identification  Scheme 

G  Polynomial 

Number  of  Parameters 


Parameter  Fixed 

Initial  Value  of  Parameters 

Q-weight i ng 

P-weight i ng 
Sampling  Interval 
Forgetting  Factor 
Initial  £c  Matrix 
Incremental  Control  Limit 


Recursive  Square  Root 
EB 

NG  1  =  3 
NF  =  2 
NH  =  1 
K 1  =1 

h  =  -1 
o 

A 1 1  set  to  zero 

kD  =  3.676  (g/s)/wt% 

=  0.919  (g/s)/wt% 

Set  to  1 

1  minute 

0.99 

1000^ 

2 . 3  g/s 
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Figure  5.13  Top  Composition  Behaviour  When  Operation  of  the 

STCWithQ-WTisInitiated 
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Figure  5.14  Variation  of  G1  Parameters  in  the  Initial 

Tuning  Period  for  STC  With  Q-WT 
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Figure  5.15  Variation  of  F  Parameters  in  the  Initial  Tuning 

Period  for  STC  With  Q-WT 
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the  system  allowed  to  reach  steady  state  and  the  flow 
increased  to  its  normal  steady  state  rate.  The  controlled 
responses  for  these  disturbances  are  shown  in  Figures  5.16 
and  5.17.  To  provide  a  quantitative  indication  of  the 
controller  performance,  the  absolute  error  values  are 
calculated  for  60  sampling  intervals  after  the  disturbance 
has  been  introduced.  These  values  are  given  in  Table  5.3. 
Also  shown  are  the  results  from  duplicate  runs. 

Table  5.3  Summary  of  Top  Composition  Absolute  Error  Values 
(Over  60  Samples)  for  STC  with  Q-weighting 


Run 

Type  of 

SAE 

No. 

F igure 

D i sturbance 

( w  t  % ) 

2-2 

5.  16 

-25%  to  Feed 

2.177 

2-3 

5.  17 

Return  to  S.S. 

2.  197 

5-4 

Dupl icate 

-25%  to  Feed 

2.347 

5-5 

Dupl i cate 

Return  to  S.S. 

1  .867 

Compar i son 

of 

the 

SAE  resu Its  with 

those  shown  in 

Table 

5.1,  shows  that 

the 

STC 

performs  better 

than 

the  PID 

control ler  for 

the 

feed 

disturbances  used. 

Whi le  the 

SAE 

achieved  by  the 

STC 

i  s 

not  much  better 

than 

the  PID  case ,  i t 

should  be  noted  that  the  PID  results  were  obtained  only 
after  considerable  time  had  been  expended  to  tune  the 
controller.  This  is  in  marked  contrast  to  the  STC  results 
obtained  using  the  first  Q-weighting  selected  (cf.  values 
given  in  Table  5.2). 
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ST  Control  of  Top  Composition  for  a  -25%  Step 
Decrease  in  Feed  Flow  Rate  With  Q-WT 


Figure  5.16 
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Figure  5.17  ST  Control  of  Top  Composition  for  a  Step 

Increase  in  Feed  Flow  Rate  to  its  Normal  Steady 
State  Value  With  Q-WT 
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5.5.2  Control  Performance  Using  the  STC  Without  O-weighting 

To  illustrate  the  effect  of  the  Q-weighting  (Q-WT),  the 
STC  was  used  to  control  the  top  composition  with  the 
Q-weighting  removed.  The  initial  controller  parameters  and 
options  are  the  same  as  those  given  in  Table  5.2  except  that 
the  Q-weighting  is  set  to  zero.  Figure  5.18  shows  the  top 
composition  response  when  the  controller  is  first  applied 
while  the  variation  of  the  parameters  as  they  adapt  during 
the  initial  identification  period  is  displayed  in  Figures 
5.19  and  5.20.  As  in  the  case  with  the  Q-weighting,  the 
parameters  have  converged  after  30  samples.  The  response  of 
the  top  composition  for  a  -25%  decrease  in  feed  flow  and  the 
subsequent  change  to  the  normal  steady  state  flow  rate  are 
shown  in  Figures  5.21  and  5.22.  The  SAE  for  these  results 
and  those  of  some  duplicate  runs  are  compared  with  the 
previous  results  obtained  using  Q-weighting  in  Table  5.4. 


Table  5.4  Comparison  of  Top  Composition  Absolute  Error 
Values  (Over  60  Samples)  for  STC  With  and 
Without  Q-weighting 


Type  of 
D i s turbance 


-25%  to  Feed 
Return  to  S.S. 

-25%  to  Feed 
Return  to  S.S. 


With 

Q-weight i ng 


SAE 

Figure  (wt%) 


5.16  2.177 

5.17  2.197 
Duplicate  2.347 
Dupl i cate  1 . 867 


Wi thout 
Q-weight i ng 


SAE 

Figure  (wt%) 


5.21  4.278 

5.22  2.693 
Duplicate  2.906 
Duplicate  3.327 
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Figure  5.18  Top  Composition  Behaviour  When  Operation  of  the 

STC  Without  Q-WT  is  Initiated 
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Variation  of  G1  Parameters  in  the  Initial 
Tuning  Period  for  STC  Without  Q-WT 


Figure  5.19 
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Figure  5.20  Variation  of  F  Parameters  in  the  Initial  Tuning 

Period  for  STC  Without  Q-WT 


■ 


90 


© 


TIME  (MIN) 


10 

•l"! 

cn 


a 
l O 
l — l 

.V  to 

:  ti 


r  u:: 

X  't 

Ti 


10 

m  .. 
rr> 


4 


6.  60 


^A/vW  '-A/ 


n - 

20.  00 


“1 

40.  00 


66.  06 


86.  00 


100 . G0 


TIME  (MEN) 


Figure  5.21 


ST  Control  of  Top  Composition  for  a  -25%  Step 
Decrease  in  Feed  Flow  Rate  Without  Q-WT 
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Figure  5.22  ST  Control  of  Top  Composition  for  a  Step 

Increase  in  Feed  Flow  Rate  to  its  Normal  Steady 
State  Value  Without  Q-WT 
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As  discussed  in  Chapter  Two,  the  Q-weighting  adds  a 
penalty  for  control  effort  so  removing  the  Q-weighting 
allows  the  manipulated  variable  to  exhibit  larger 
deviations.  Comparison  of  the  variations  in  reflux  flow 
rate  shown  in  Figure  5.21  and  5.16  shows  that  this  is  indeed 
the  case.  It  can  be  seen  that  without  Q-weighting,  the 
magnitude  of  the  reflux  flow  variation  is  twice  the 
variation  that  resulted  when  Q-weighting  was  used.  The 
larger  variation  in  reflux  flow  due  to  the  removal  of  the 
Q-weighting  causes  the  control  performance  to  deteriorate  as 
shown  by  the  SAE  values  for  STC  with  and  without  Q-weighting 
given  in  Table  5.4. 

Comparing  the  initial  tuning  in  period  for  the  cases 
with  and  without  Q-weighting,  it  is  seen  that  the  added 
penalty  in  control  effort  helps  to  reduce  the  variation  in 
reflux  flow.  In  addition,  without  the  Q-weighting,  the 
controller  demands  a  maximum  reflux  flow  of  over  16  g/s  but 
using  the  Q-weighting,  the  maximum  is  reduced  to  11.5  g/s. 
This  reduction  in  the  maximum  control  output  leads  to  less 
overshoot  in  the  overhead  composition  response. 

These  results  clearly  show  that  the  use  of  the 
Q-weighting  reduces  the  amount  of  variation  in  the  reflux 
flow  rate.  This  will  lead  to  an  improvement  in  the 
controller  performance  in  terms  of  the  SAE  criteria. 
Futhermore,  the  use  of  Q-weighting  does  not  slow  down  the 
rate  at  which  the  parameters  converge  during  the  initial 
tuning  in  period  and  it  gives  the  added  advantage  of 
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reducing  the  maximum  reflux  flow  required  during  the  initial 
transition  period. 

5.5.3  Effect  of  Sampling  Time 

To  show  the  effect  of  sampling  time,  experimental  runs 
were  made  using  times  of  two  and  three  minutes.  The 
experimental  responses  for  a  -25%  step  decrease  in  feed  flow 
rate,  for  the  two  and  three  minute  rates  are  shown  in 
Figures  5.23  and  5.24.  With  the  exception  of  the  change  in 
sampling  time  and  the  absence  of  the  Q-weighting,  the 
initial  controller  parameters  and  options  are  the  same  as 
those  shown  in  Table  5.2.  No  Q-weighting  was  employed  in 
order  to  demonstrate  that  increasing  the  sampling  time  has 
the  same  effect  as  adding  the  Q-weighting;  that  is  in  terms 
of  reducing  the  oscillation  in  the  manipulated  var i able[ 24 ] . 

Comparing  the  overhead  composition  response  for  the 
different  sampling  times,  it  is  seen  that  by  increasing  the 
time  to  two  minutes,  the  magnitude  of  the  variation  in  the 
reflux  flow  is  greatly  reduced.  As  can  be  seen  from  the 
absolute  error  values  in  Table  5.5,  the  performance  is  not 
much  better  than  that  achieved  using  the  one  minute  sampling 
rate  because  at  the  two  minute  sampling  time,  the  controller 
responds  more  slowly  to  an  unknown  disturbance  since  the 
controller  is  only  acting  in  the  feedback  mode.  By  further 
increasing  the  sampling  time  to  three  minutes,  the  magnitude 
of  the  oscillation  in  the  reflux  flow  starts  to  increase. 

The  controlled  response  to  the  disturbance  is  not  as 
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Figure  5.23  ST  Control  of  Top  Composition  to  a  -25%  Step 

Decrease  in  Feed  Flow  Rate  Without  Q-Wt:  Two 
Minute  Sampling  Time 
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Fiqure  5.24  ST  Control  of  Top  Composition  to  a  -25%  Step 

Decrease  in  Feed  Flow  Rate  Without  Q-Wt:  Three 
Minute  Sampling  Time 
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satisfactory  because  of  the  larger  sampling  time. 


Table  5.5  Comparison  of  Top  Composition  Absolute  Error 

Values  Calculated  for  a  Period  of  One  Hour  After 
the  Feed  Disturbance  for  Different  Sampling  Times 


Sampling  Time  SAE 
Figure  (minutes)  (wt%) 


5.16  1  4.278 

5.23  2  3.627 

5.24  3  6.700 


The  results  show  that  while  it  is  possible  to  use 
larger  sampling  times  to  reduce  the  variation  in  the  reflux 
flow,  the  larger  sampling  time  tends  to  slow  down  the 
controller  response  for  a  disturbance  that  has  to  be 
detected  by  feedback  of  the  output  variable.  For  a  feedback 
strategy,  it  would  be  better  to  sample  faster  and  increase 
the  weighting  placed  on  the  control  effort.  However,  for  a 
controller  with  feedforward  action,  it  should  be  possible  to 
achieve  the  same  control  performance  with  a  larger  sampling 
time. 

5.5.4  Effect  on  Control  Performance  when  Identification  is 
Stopped 

From  Equation  (2.34),  the  control  law  can  be  expressed 
as : 

Fy  +  Gu  -  HQu  +  Hw^  +  Dv  =0  (5.1) 

t+k  t  t  t  t 

and  without  feedforward  action,  the  control  law  reduces  to: 
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Fy  +  Gu  -  HQu  +  Hw  =  0  (5.2) 

t+k  t  t  t 

At  steady  state,  the  outputs  and  inputs  remain  constant 
and  Qu  tends  to  zero  (cf.  Chapter  2).  Equation  (5.2) 
becomes : 

XFiyt+k  +  EGiut  '  EHist  =  0  (5-3) 

If  the  identification  scheme  is  stopped,  the  sum  of  the 

coefficients  for  each  of  the  polynomials  becomes  constant. 

When  a  disturbance  is  introduced,  it  is  clear  that  the 

control  output  will  reach  a  new  steady  state  value.  Since 

the  sum  of  these  coefficents  remains  the  same,  to  satisfy 

Equation  (5.3),  only  the  process  output  value  can  change 

hence  resulting  in  an  offset.  The  magnitude  of  the  offset 

will  depend  on  the  G  and  F  parameters.  If  the  sum  of  the  G 

parameters  is  small  compared  to  the  sum  of  the  F  parameters, 

only  a  small  change  in  y  is  needed,  so  a  smaller  offset 

t  k 

occurs . 

To  demonstrate  the  effect  discussed,  a  run  was  made 
where  the  identification  was  stopped  after  the  STC  had 
brought  the  top  composition  to  its  steady  state  value.  The 
initial  controller  parameters  and  options  are  the  same  as 
those  given  in  Table  5.2  with  the  exception  that  Q-weighting 
is  not  used.  The  response  of  the  top  composition  to  a  -25% 
step  change  from  the  steady  state  feed  flow  rate  is  recorded 
in  Figure  5.25.  Although  the  offset  is  small,  it  is  clear 
that  there  is  an  offset  of  about  0.05%.  By  resuming  the 
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Figure  5.25  ST  Control  of  Top  Composition  to  a  -25%  Step 

Change  in  Feed  Flow  Rate  Without  Q-WT  and  the 
Identification  Stopped 
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indenti f ication,  the  parameters  will  readjust  to  remove  the 
offset  as  can  be  seen  from  the  response  shown  in  Figure 
5.26. 

5.6  BOTTOM  COMPOSITION  CONTROL  WITH  STC 

The  self-tuning  controller  was  used  to  control  the 
bottom  composition,  despite  changes  in  feed  flow  rate,  by 
manipulating  the  steam  flow  rate.  The  reflux  flow  was  fixed 
at  its  normal  steady  state  value  of  10.2  g/s.  The  initial 
self-tuning  controller  parameters  and  options  are  given  in 
Table  5.6.  The  control  behavior  for  bottom  composition  that 
resulted  for  the  four  different  feed  flow  rate  disturbances 
used  in  this  study  are  shown  in  Figures  5.27  to  5.30.  The 
results  given  are  for  STC  runs  that  already  have  'good' 
parameter  estimates.  Some  re-tuning  of  the  Q-weighting  was 
necessary  with  the  resulting  values  as  noted.  Table  5.7 
gives  a  summary  of  the  absolute  error  values  for  bottom 
composition  control  for  a  period  of  1  1/2  hour  after  the 
feed  rate  was  changed  for  both  the  STC  and  the  previous  PID 
results  (cf.  Table  5.1).  It  can  be  seen  from  these  values 
that  the  STC  performs  better  than  PID  for  the  +25%  increase 
in  feed  flow  rate  away  from  the  steady  state  value  and  for 
the  decrease  back  to  steady  state  feed  flow,  while  PID  is 
better  than  the  STC  for  the  other  two  disturbances.  On  the 
basis  of  these  results,  it  can  be  seen  that  both  controllers 
performed  equally  well.  However,  it  should  be  noted  that 
considerable  more  time  was  needed  to  tune  the  PID 
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Figure  5.26 


Top  Composition  Response  After  Resuming 
Identification  in  the  STC 


Table  5.6  Controller  Parameters  and  Options  for  STC 
Bottom  Composition  Control 


Identification  Scheme 

G  Polynomial 

Number  of  Parameters 


Parameter  Fixed 

Initial  Value  of  Parameters 

Q-weight i ng 

P -weight i ng 
Sampl i ng  Interva 1 
Forgetting  Factor 
Initial  £c  Matrix 
Incremental  Control  Limit 


Recursive  Square  Root 
EB 

NG1  =  4 
NF  =  3 

ND  =  5  (FF  added) 

NH  =  1 
K 1  =1 

K3  =  1  (FF  added) 


All  set  to  zero 

K  =  1.145  ( g/s ) /wt% 
k.  =  0.226  ( g / s ) / w t % 

Set  to  1 

3  minute 

0.99 

1000^ 

2.25  g/s 
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Figure  5.27  ST  Control  of  Bottom  Composition  for  a  -25% 

Step  Change  in  Feed  Flow  Rate  With  Q-WT 
( K  p  -  1.15  and  K  j  -  0.23) 
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Figure  5.28  ST  Control  of  Bottom  Composition  for  a  Step 

Increase  in  Feed  Flow  Rate  to  its  Normal  Steady 
State  Value  With  Q-WT  (Kp  -  0.90  and  -  0.47) 
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g  e  5.29  ST  Control  of  Bottom  Composition  for  a  +25% 
Step  Change  in  Feed  Flow  Rate  With  Q-WT 
( K  -0.90  and  k  -0.47) 
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Figure  5.30  ST  Control  of  Bottom  Composition  for  a  Step 

Decrease  in  Feed  Flow  Rate  to  its  Normal  Steady 
State  Value  With  Q-WT  (Kp  -  0.90  and  Kj  -  0.47) 
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compensator  than  the  Q-weighting  constants.  Also,  the  final 
Q-weighting  parameters  are  closer  to  the  original  estimates. 
In  addition,  the  PID  constants  for  the  different  disturbances 
differ  by  a  considerable  amount  in  comparison  to  the 
variations  in  the  Q-weighting  parameters  for  the  different 
disturbances.  This  can  be  seen  by  comparing  the  difference 
in  the  two  sets  of  PID  constants  for  bottom  composition 
control  given  in  Table  5.1  with  the  difference  in  the 
Q-weighting  parameters  for  the  responses  shown  in  Figure 
5.27  to  5.30.  This  is  because  the  STC  is  able  to  adapt  to 
the  nonlinear  dynamic  behavior  of  the  column. 

Table  5.7  Comparison  of  Bottom  Composition  Absolute  Error 
Values  (over  1  1/2  hour)  for  PID  and  STC 

PID  Control  ST  Control 


Type  of 

Di sturbance 

Figure 

SAE 

(wt%) 

Figure 

SAE 

(wt%) 

-25%  to  Feed 

5.7 

18.43 

5.27 

19.70 

Return  to  S.S. 

5.8 

14.00 

5.28 

14.90 

+25%  to  Feed 

5.5 

20.72 

5.29 

18.76 

Return  to  S.S. 

5.6 

26.04 

5.30 

20.51 

Additional  evidence  that  the  STC  handles  nonlinearity 
better  than  the  PID  compensator  is  seen  by  comparing  the 
performances  of  the  two  controllers  for  the  +25%  step 
increase  in  feed  flow  and  the  subsequent  step  decrease  to 
the  normal  feed  rate  (cf.  Figures  5. 6-5. 7  and  5.29-5.30). 
From  the  SAE  values  given  in  Table  5.7,  it  can  be  seen  that 
the  STC,  compared  to  the  PID  compensator,  provided  similar 
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control  performance  for  theses  different  step  changes  in 
feed  flow. 

The  effect  of  the  addition  of  feedforward  (FF)  control 
action  to  the  STC  was  studied  using  the  same  controller 
parameters  and  options  given  in  Table  5.6  except  that  the 
values  of  the  G1  and  F  parameters  are  initially  set  equal  to 
the  final  values  obtained  when  using  the  STC  without  feed¬ 
forward  action. 

The  bottom  composition  responses  for  the  same  four 
different  feed  flow  rate  disturbances  are  given  in  Figures 
5.31  to  5.34  with  the  absolute  error  values  for  these 
responses  and  those  of  the  previous  PID  and  STC  tests 
summarized  in  Table  5.8. 

Table  5.8  Comparison  of  Bottom  Composition  Absolute  Error 
Values  (over  1  1/2  hour)  for  PID,  STC  and  STC+FF 


Type  of 

Di s turbance 

PID  Control 

ST  Control 

STC 

+  FF 

F igure 

SAE 

(wt%) 

F igure 

SAE 

( w  t  %  ) 

F igure 

SAE 

( w  t  %  ) 

-25%  to  Feed 

5.7 

18.43 

5.27 

19.70 

5.31 

26.27 

Return  to  S.S. 

5.8 

14.00 

5.28 

14.90 

5.32 

9.93 

+25%  to  Feed 

5.5 

20.72 

5.29 

18.76 

5.33 

9.  17 

Return  to  S.S. 

5.6 

26.04 

5.30 

20.51 

5.34 

8.41 

The  results  given 

are  for 

STC  runs  that 

a  1  ready 

have 

'good'  parameter  estimates.  Except  for  the  -25%  step 
decrease  in  feed  flow  rate  from  its  steady  state  value,  the 
STC  with  feedforward  control  action  improves  the  bottom 


XB  -  CUTX3  ST  -  LG/ S3 

.  0J  4.00  7.0 3  10.00  11.00  13.00  15.00  17.00 


108 


0.00  30.00  60.  00  30.  00  120.00  130.00 

TIME  (MIN) 


TIME  (MIN) 


Figure  5.31  STC  +  FF  Control  of  Bottom  Composition  for  a 

-25%  Step  Change  in  Feed  Flow  Rate  With  Q-WT 
(Kp  -  0.90  and  -  0.47) 
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STC  +  FF  Control  of  Bottom  Composition  for  a 
Step  Increase  in  Feed  Flow  Rate  to  its  Normal 
Steady  State  Value  With  Q-WT  (kD  -  0.90  and 
kj  -  0.47) 


Figure  5.32 
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( Kp  -  0.90  and  K  s  -  0.47) 
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Figure  5.34 


composition  control  significantly.  The  error  values  using 
feedforward  are  almost  one-half  of  those  values  obtained 
when  feedforward  control  action  was  not  used.  The 
self-tuning  controller  did  not  perform  better  than  the  PID 
controller  for  the  -25%  step  decrease  from  the  steady  state 
feed  flow  rate  because  the  feedforward  terms  have  not  had 
time  to  be  identified. 


5.7  SIMULTANEOUS  CONTROL  OF  TERMINAL  COMPOSITION 
Simultaneous  control  of  overhead  and  bottom 
compositions  using  the  self-tuning  controller  was  carried 
out  for  each  of  the  different  feed  disturbances.  The 
terminal  compositions  were  controlled  using  two  separate  STC 
(multi  loop  conf igurat ion , ML-STC )  and  two  STC  that  were 
linked  by  feedforward  terms  to  cancel  the  effect  of 
interaction  (multivariable  conf igurat ion , MV-STC ) . 

5.7.1  Multi  loop  Self-Tuning  Control  of  Terminal  Composition 
The  original  controller  parameters  and  options  for  the 
two  separate  self-tuning  controllers  are  given  in  Table  5.9. 
The  intial  parameters  values  are  the  final  estimates  from 
the  single  composition  control  tests. 

After  an  intial  identification  period,  the  feed 
disturbance  is  introduced.  It  was  found  that  the  top  loop 
Q-weighting  had  to  be  de-tuned  to  prevent  oscillatory 
behavior.  The  composition  behavior,  for  the  de-tuned 
Q-weighting  test  for  the  -25%  step  decrease  from  the  steady 
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Table  5.9  Controller  Parameters  and  Options  for  Multi  loop 
and  Multivariable  STC 


Recursive  Square  root  (top  and  bottom) 
EB  (top  and  bottom) 


Ident i f i ca t i on 
Scheme 

G  Polynomial 

Number  of 
Parameters 


Top 


NG 1 

3 

NF 

— 

2 

ND 

z 

0 

NH 

— 

1 

NG2 

— 

2 

(MV-STC) 

K  1 

z 

1 

K2 

z 

0 

(MV-STC) 

K3 

zz 

0 

Bot  tom 
NG1  =  4 
NF  =  3 

ND  =  5  (FF  added) 
NH  =  1 

NG2  =  3  (MV-STC) 

K 1  =  2 

K2  =  0  (MV-STC) 

K3  =  1  (FF  added) 


Parameter  fixed 


h  =  -1  (top  and  bottom) 
o 


Initial  Value  of 
Parameters 


Top 

G 1  F 

0.2090  0.7382 

-  .0409  0.1161 

-  .0834 


Bot  tom 

G 1  F  D 

0.1586  1.0340  0.03257 

0.0548  0.1695  0.04689 

-.0997  -.6200  0.00937 

-.0763  -.03920 

- .03590 


Q-weighting  : 

P-weighting  : 

Sampling  Interval: 


Top 

kn  =  3.676 
KT  =  0.919 

Set  to  1  (top 

3  mi nute  ( top 


Bot  tom 
1  .  145 
0.226 

and  bottom) 
and  bottom) 


( g/s ) /wt% 
( g/s ) /wt% 


Forgetting  Factor:  0.99  (top  and  bottom) 


Initial  : 

£c  Matrix  :  0.011  (top  and  bottom) 

Incremental 

Control  Limit  :  Top  Bottom 

2.30  g/s  2.25  g/s 
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state  feed  flow  rate  and  the  subsequent  return  to  the  normal 
value  are  shown  in  Figures  5.35  and  5.36.  The  control 
performance  that  was  obtained  with  feedforward  control 
action  added  to  the  bottom  STC  is  shown  in  Figures  5.37  and 
5.38.  A  summary  of  the  absolute  error  values  over  a  period 
of  30  samples  after  the  feed  flow  rate  disturbance  was 
introduced  is  given  in  Table  5.10. 

Table  5.10  Summary  of  Top  and  Bottom  Composition  Absolute 

Error  Values  for  Multi  loop  STC  With  and  Without 
Feedforward  (FF)  added 


Type  Of  Disturbance 
-25%  to  Feed  Return  to  S.S.  +25%  to  Feed 


Type  of 
Control 

F igure 

SAE 

<wt%) 

Figure 

SAE 

( w  t  % ) 

F igure 

SAE 

(wt%) 

ML  -  STC 

top 

bot 

5.35 

5.35 

2.32 

23.75 

5.36 

5.36 

2.  17 
19.71 

ML  -  STC 
+  FF 

top 

bot 

5.37 

5.37 

3.74 

14.59 

5.38 

5.38 

2.27 

8.97 

The  SAE  values  show  that  the  control  performance  of  the 
multi  loop  self-tuning  controller  was  comparable  to  that  of 
the  well  tuned  PID  compensator.  It  should  be  noted  that  as 
with  bottom  composition  control,  implementation  of  the  STC 
required  much  less  tuning  than  was  the  case  for  the  PID 
regulator.  In  fact,  the  Q-weighting  for  the  bottom 
composition  STC  was  not  re- tuned,  from  the  single  loop 
values,  for  use  with  multi  loop  control.  In  addition, 
proportional,  integral  and  derivative  actions  are  involved 
in  the  PID  compensator  while  only  proportional  and  integral 
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ML-ST  Control  of  Terminal  Compositions  for  a 
Step  Increase  in  Feed  Flow  Rate  to  its  Normal 
Steady  State  Value  With  Q-WT  (top:  Kp  -  2.7 
Kj  -  .45bot:  kp  -  .90  K(  -  .47) 


Figure  5.36 
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Figure  5.37  ML-STC  +  FF  (bottom)  Control  of  Terminal 

Compositions  for  a  Step  Increase  in  Feed  Flow 
Rate  to  its  Normal  Steady  State  Value  With  Q-WT 
(top:  Kp  -  2.7  kj  -  .45  bot :  kp  -  .90  k(  -  .47) 
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Figure  5.38  ML-STC  +  FF  (bottom)  Control  of  Terminal 

Compositions  for  a  +25%  Step  Change  in  Feed 
Flow  Rate  With  Q-WT  (top:  Kp  -  2.7  Kj  -  .45 
bot :  Kp  -  . 90  Kj  -  . 47 ) 
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actions  are  used  with  the  Q-weighting  of  the  self-tuning 
controller.  The  STC  response  should  improve  if  derivative 
action  is  added. 

The  addition  of  feedforward  action  for  the  bottom 
self-tuning  controller  has  significantly  improved  the 
control  performance  as  can  be  seen  by  examining  the  SAE 
values  in  Table  5.10.  The  bottom  composition  SAE  values  are 
greatly  reduced  while  the  top  composition  SAE  values  for  the 
STC  are  comparable  to  the  values  that  result  with  PID 
control . 

5.7.2  Multivariable  ST -Control  of  Terminal  Composition 

A  multivariable  STC  was  also  used  to  control  the 
terminal  compositions.  The  original  parameters  and  options 
used  by  the  multivariable  STC  are  the  same  as  those  in  Table 
5.8  except  that  the  feedforward  linkage  terms  (G2)  are 
identified.  The  composition  responses  for  the  25%  step 
decrease  in  feed  flow  rate  and  the  subsequent  increase  in 
feed  flow  rate  to  return  the  flow  to  the  normal  steady  state 
value  are  given  in  Figures  5.39  and  5.40. 

The  absolute  error  values  for  these  responses  and  those 
of  the  multi  loop  STC  and  PID  control  are  given  in  Table 
5.11.  These  error  values  show  that  the  multivariable  STC 
performed  as  well  as  the  PID  compensator  for  bottom 
composition  control,  but  provided  significantly  better 
control  of  top  composition  control  than  did  PID  control. 

The  top  composition  error  values  with  multivariable  STC  are 
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Figure  5.39 
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122 


about  one-half  the  values  that  result  using  PID  control. 
The  improvement  is  due  to  the  addition  of  the  feedforward 
terms  which  link  the  two  controllers  and  serve  to  reduce 
interaction  hence  improving  control. 

Table  5.11  Comparison  of  Top  and  Bottom  Composition 

Absolute  Error  Values  for  Multi  loop  PID, 

Multi  loop  STC  (With  and  Without  F F ) ,  and 
Multivariable  STC 


T  y 

p  e  0 

f  D  i  s 

t  u  r 

banc 

e 

-25%  to 

Feed 

Return  to 

S.S. 

+25%  to 

Feed 

Type  of 
Control 

F igure 

SAE 
(wt% ) 

Figure 

SAE 
( wt% ) 

Figure 

SAE 

(wt% 

ML  -  PID 

top 

bot 

5.11 

5.11 

2.13 

21.69 

5.  12 

5.  12 

2.30 

22.96 

5.9 

5.9 

2.75 

27.26 

ML  -  STC 

top 

bot 

5.35 

5.35 

2.32 

23.75 

5.36 

5.36 

2.17 

19.71 

ML  -  STC 
+  FF 

top 

bot 

5.37 

5.37 

3.74 

14.59 

5.38 

5.38 

2.27 

8.97 

MV  -  STC 

top 

bot 

5.39 

5.39 

1.38 

27.39 

5.40 

5.40 

1  .07 
17.28 

5.8  DISCUSSION  OF  RESULTS 

From  the  results  given,  it  has  been  shown  that  the 
feedback  self-tuning  controllers  can  control  the  terminal 
compositions  of  the  distillation  column  as  well  as  or  better 
than  well  tuned  PID  compensators.  With  the  multivariable 
STC,  the  top  composition  controller  is  able  to  react 
immediately  to  steam  changes  hence  improving  the  top 
composition  control  compared  to  the  multi  loop  PID  case. 
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The  ability  of  the  feedback  STC  to  provide  a  control 
performance  that  is  comparable  to  a  PID  compensator  is 
remarkable  considering  the  type  of  disturbance  used  to 
compare  the  controllers.  The  use  of  step  changes  in  feed 
flow  rate  as  disturbances  is  a  very  severe  test.  In  the 
derivation  of  the  STC,  the  disturbances  have  been  assumed  to 
have  zero  mean.  When  a  step  change  in  the  feed  flow  rate  is 
introduced,  a  new  steady  state  control  output  is  required, 
as  calculated  from  the  control  law.  Therefore,  the  feed 
flow  step  changes  are  not  zero  mean  disturbances.  To 
further  appreciate  this  situation,  consider  what  happens  to 
the  STC  when  the  feed  rate  is  increased.  As  the  additional 
flow  of  material  moves  through  the  column,  the  bottom 
methanol  composition  starts  to  increase  and  the  STC  compares 
the  composition  with  the  predicted  value  and  finds  there  is 
an  error.  The  controller  will  then  increase  the  steam  flow 
rate  and  in  addition  change  the  parameters  so  that  they  will 
give  a  prediction  consistent  with  the  higher  recorded  value. 
As  more  of  the  increased  feed  reaches  the  reboiler,  the 
composition  becomes  even  higher,  the  steam  flow  rate  is 
again  increased  and  the  parameters  further  changed  to  give 
the  correct  output.  In  effect,  there  is  a  transition  period 
when  the  parameters  change  and  the  predicted  outputno  long¬ 
er  matches  the  actual  output.  Figures  5.41  and  5.42  show 
how  the  bottom  STC  parameters  change  when  a  -25%  step  change 
in  feed  flow  rate  occurs  (run  shown  in  Fig  5.27).  It  is 
seen  that  there  is  a  significant  change  in  the  parameter 
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Figure  5.41  Variation  of  G1  Parameters  in  the  Transition  Period  of  a  -25%  Step 

Change  in  Feed  Flow  Rate  for  the  Bottom  5TC  with  Q-WT 
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Figure  5.42  Variation  of  F  Parameters  in  the  Transition 

Period  of  a  -25%  Step  Change,  in  Feed  Flow  Rate 
for  the  Bottom  STC  with  Q-WT 
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estimates.  It  is  due  to  the  robustness  of  the  STC  that  using 
only  feedback  information,  it  is  able  to  match  the  control 
performance  of  the  well  tuned  PID.  To  properly  handle  step 
changes  in  feed  flow  rate,  the  feedforward  option  should  be 
used.  The  results  of  the  self-tuning  controller  with 
feedforward  added  does  indeed  show  considerable  improvement 
over  the  PID  case. 

It  should  also  be  noted  that  the  PID  control  behavior 
shown  has  been  obtained  after  considerable  time  was  spent 
fine  tuning  the  compensator.  On  the  other  hand,  the  amount 
of  time  needed  to  tune  the  Q-weighting  was  minimal.  The 
original  estimates  obtained  from  Halman's  technique  [75] 
gave  good  control.  With  more  tuning  of  the  Q-weighting,  it 
should  be  possible  to  improve  the  control  performace. 
Additional  improvement  in  the  controller  performance  should 
also  be  realized  if  derivative  action  is  added  into  the 
Q-weighting . 
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6.  CONCLUSIONS  AND  RECOMENDAT IONS 


This  work  was  concerned  primarily  with  the  experimental 
evaluation  of  the  self-tuning  controller.  However,  before 
the  evaluation  could  be  carried  out,  the  pairing  of 
manipulated  and  controlled  variables  had  to  be  decided.  The 
main  concern  was  whether  top  composition  should  be 
controlled  by  manipulating  reflux  flow  (energy  balance 
control)  or  distillate  flow  (material  balance  control).  It 
was  found  that  compared  to  energy  balance  control,  the 
material  balance  control  strategy  is  not  as  sensitive  to 
heat  disturbance  provided  the  external  reflux  ratio  is 
large.  However  just  the  reverse  case  is  true  if  the 
external  reflux  ratio  is  small.  Analysis  of  the  sensitivity 
of  the  vapor- to- 1 iquid  ratio  to  changes  in  the  vapor  rates 
for  the  two  strategies  provides  valuable  information  on 
which  of  the  two  control  schemes  should  provide  better 
control  performance  for  heat  disturbances  entering  the 
column.  For  the  pilot  plant  distillation  column  used  in 
this  work,  analysis  of  the  vapor- to- 1 iquid  ratio  indicated 
that  energy  balance  control  was  superior  and  this  was 
substantiated  by  the  results  from  experimental  tests. 

To  evaluate  the  STC,  a  comparison  was  made  between  the 
performance  of  the  STC  and  a  well  tuned  PID  compensator  in 
controlling  the  terminal  compositions  of  the  pilot  plant 
column  when  it  was  subjected  to  step  changes  in  the  feed 
rate.  The  evaluation  was  carried  out  for:  top  composition 
control,  bottom  composition  control  and  dual  composition 
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control.  In  all  cases,  it  was  found  that  the  STC  is  superior 
to  the  PID  compensator.  Specifically,  it  has  been  shown 
that : 

(1)  The  STC  acting  only  in  the  feedback  mode  using  initial 
Q-weighting  estimates  from  Halman's  [75]  technique 
performed  as  well  as  a  PID  compensator  that  had  been 
very  well  tuned  at  a  considerable  expenditure  of  time. 

(2)  Addition  of  feedforward  control  action  on  the  bottom 
composition  loop  significantly  improved  the  controller 
per  formance . 

(3)  Use  of  the  multivariable  STC  allowed  the  top 
composition  controller  to  react  immediately  to  changes 
in  heat  input  hence  improve  top  composition  control 
per formance . 

In  addition,  the  STC  was  found  to  have  the  following 
proper  ties: 

(1)  Starting  with  zero  estimates,  the  STC  can,  for  a  fast 
system  (top  composition  loop),  converge  to  reasonable 
estimates  after  only  30  sampling  intervals.  It  gives 
superior  control  compared  to  PID  control  results  when 
the  STC  is  used  with  these  estimates. 

(2)  Fixing  the  hQ  parameter  value  to  -1  does  not  create  any 
problems  with  parameter  convergence. 

(3)  The  addition  of  Q-weighting  reduces  the  variation  in 
the  manipulated  variables  and  as  a  result  the  control 
of  the  output  variable  is  improved  considerably. 

The  controller  parameters  are  able  to  readjust  to  some 
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degree  for  changes  in  the  Q-weighting  so  that  the 
output  response  is  not  as  sensitive  to  changes  in  the 
Q-weighting  constants  as  in  the  case  of  the  PID 
compensator.  Nevertheless,  the  Q-weighting  can  still 
be  fine  tuned  the  same  way  as  a  normal  PI  compensator 
is  tuned  to  obtain  even  improved  control  responses. 

The  relationship  between  K  and  K.  and  how  they  affect 

pi  y 

the  output  response  remains  the  same. 

(5)  If  the  identification  phase  of  the  STC  is  stopped,  an 
unknown  disturbance  will  cause  an  offset  to  occur.  The 
magnitude  of  the  offset  depends  upon  the  sum  of  the  F 
and  G  parameters. 

(6)  To  an  extent,  a  large  sampling  time  has  the  same  effect 
as  adding  more  Q-weighting;  i.e.,  it  smooths  out  the 
variation  in  the  manipulated  variable.  However,  if  the 
STC  is  acting  only  in  the  feedback  mode,  the  controller 
response  will  suffer  because  it  reacts  more  slowly  to 
unknown  disturbances  due  to  the  larger  sampling  time. 

In  summary,  the  experimental  evaluation  of  the 
self-tuning  controller  as  modified  by  Morris  [24]  has  shown 
that  the  controller  is  superior  to  the  PID  compensator.  If 
tuning  is  necessary,  the  controller  can  be  fine  tuned  the 
same  way  a  PID  regulator  is  tuned.  The  controller  is 
extremely  robust  as  it  is  able  to  handle  non-zero  mean 
disturbances  even  though  the  control  law  was  derived  based 
on  the  assumption  that  the  disturbance  is  a  white  noiise 
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From  the  results  of  this  work,  the  following 

recommendations  are  given  for  possible  future  studies: 

(1)  With  the  present  column  operating  conditions,  energy 
balance  control  has  been  shown  to  be  superior.  It  is 
recommended  that  test  be  performed  with  the  operating 
conditions  altered  so  that  the  external  reflux  ratio  is 
higher.  This  would  indicate  whether  material  balance 
control  will  indeed  be  less  sensitive  to  heat 
disturbances  under  the  new  operating  conditions. 

(2)  While  performing  the  experimental  work  for  this  study, 
observations  on  the  multipoint  temperature  chart 
recorder  indicated  that  the  liquid  temperature  on  tray 
1  gave  a  good  indication  of  changes  in  bottom 
composition.  It  is  suggested  that  bottom  composition 
be  controlled  by  maintaining  tray  1  temperature 
constant  using  the  self-tuning  controller.  For  a 
binary  system,  providing  the  pressure  is  constant, 
fixing  the  temperature  implies  that  the  composition 
will  remain  constant.  To  account  for  any  variation  in 
pressure,  the  result  of  the  GC  analysis  can  be  used  to 
reset  the  tray  1  temperature  set  point. 

(3)  The  self-tuning  controller  is  designed  primarily  to 
handle  stochastic  disturbances  while  the  evaluation  has 
been  carried  out  using  deterministic  disturbances.  A 
further  evaluation  should  be  made  on  how  the 
performance  of  the  STC  compares  with  the  PID 
compensator  when  random  disturbances  are  used;  e.g., 
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add  noise  to  the  composition  measurements. 

(4)  When  an  unknown  deterministic  disturbance  is 
encountered,  the  STC  forces  the  parameters  to  change 
rapidly  to  account  for  the  disturbance.  An  improved 
response  might  be  achieved  by  freezing  the  parameters 
once  reasonable  estimates  have  been  obtained.  To 
prevent  the  offset  that  occurs  when  i dent i f i caton  is 
stopped,  an  integrator  can  be  added  to  the  controller. 
The  control  action  would  be  the  normal  STC  algorithm 
output  plus  an  added  action  based  on  the  error  from  set 
point. 

(5)  In  the  present  implementation  of  the  STC  algorithm,  the 
feedforward  control  action  option  only  allows  the  use 
of  lead  terms.  As  lag  action  is  normally  required  in 
feedforward  compensation,  it  should  be  added  into  the 
algor i thm. 

(6)  This  study  considered  the  performance  of  the 
self-tuning  control ler  only  for  load  disturbances. 
Additional  testing  should  be  performed  to  study  the 
performance  of  the  control ler  for  set  point  changes. 

(7)  The  multivariable  approach  given  by  Morris  [24]  which 
is  used  in  this  study  differs  from  that  proposed  by 
Keviczky  [19]  and  Borrison  [19]  in  that  the 
identification  of  the  system  is  performed  on  a  loop  by 
loop  basis.  Both  Keviczky  and  Borrison  suggested 
identifying  the  multivariable  model  directly.  A  direct 
comparison  by  simulation  or  experimentally  of  these 
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different  multivariable  self-tuning  schemes  should  be 
undertaken . 

(8)  Rather  than  starting  with  zero  initial  estimates,  it  is 
suggested  that  the  STC  be  allowed  to  identify  para¬ 
meter  estimates  while  the  process  is  operating  under 
PID  control.  These  estimates  should  provide  a 

smooth  transition  during  the  initial  switch  to  the 
self-tuning  controller. 

(9)  The  problem  of  the  covariance  matrix  bursting  due  to 
inadequate  excitation  should  be  investigated.  The 
different  methods  that  have  been  proposed  to  overcome 
the  problem  should  be  tested;  eg.  use  of  recursive 
learning  identification  as  proposed  by  Morris  [24]. 
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NOMENCLATURE 


a.  Abbreviations: 

DACS  -  Data  Acquisition  Control  and  Simulation 


DDC 

-  Direct  Digital  Control 

FF 

-  FeedForward 

GC 

-  Gas  Chromatograph 

HP 

-  Hewlett  Packard 

ML 

-  Multi  loop 

MV 

-  Multivariable 

PID 

-  Proportional -Integral -Derivative 

Q-WT 

-  Q -weigh ting 

SAE 

-  Sum  of  Absolute  Error 

SISO 

-  Single  Input  Single  Output 

STR 

-  Self-Tuning  Regulator 

ST  C 

-  Self-Tuning  Controller 

S.S. 

-  Steady  State 

Variables: 


D 

Distillate  flow  rate 

FE  - 

Feed  flow  rate 

d 

Cost  function  for  SISO  case 

J 

Cost  function  for  multivariable  case 

I 

Ident i ty  matrix 

k 

System  time  delay  in  integer  sample  intervals. 

kA  - 

Time  delay  due  to  other  output  for  a  multivariable 
system 
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D 

kD  -  Time  delay  due  to  other  control  input  for  a 
multivariable  system 

k*-  -  Time  delay  associated  with  load  disturbance 

kp  -  Proportional  constant  in  PID  compensator 

k.  -  Integral  gain  in  PID  compensator 

k^  -  Derivative  gain  in  PID  compensator 

K  -  Kalman  gain  vector  used  in  parameter  estimation 
algor i thm 

K1  -  Use  in  STC  program  to  designate  time  delay  of 
system 

K2  -  Use  in  STC  program  to  designate  time  delay  due  to 

other  control  input  for  a  multivariable  system: 

kL  -  kB 
•  •  •  • 

i  j  n 

K3  -  Use  in  STC  program  to  disignate  time  delay 
associated  with  load  disturbance:  kB  -  k 

L  -  Liquid  rate  inside  the  column 

M  -  Number  of  loops  in  multivariable  system 


n 

-  Sys 

tern 

order 

ND 

-  No. 

of 

D  parameters 

i  n 

control  1 er 

NG 1 

-  No. 

of 

G1  parameters 

i  n 

control ler 

NG2 

-  No . 

of 

G2  parameters 

i  n 

control ler 

NF 

-  No . 

of 

F  parameters 

i  n 

control ler 

NH 

-  No . 

of 

H  parameters 

i  n 

control ler 

pc 

-  Covari 

ance  Matrix 

RE 

-  Ref 

lux 

flow  rate 

ST 

-  Steam 

flow  rate 

t  -  Time  measured  in  integral  sample  intervals 

u  -  Control  input 

V  -  Vapor  rate  inside  the  column 

v  -  Load  disturbance 
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v  -  Modified  load  disturbance  defined  as  zk"kLv 
w  -  Set  point 

w  -  Modified  set  point  defined  as  Hw 
X  -  Vector  containing  measured  data 
XB  -  Bottom  Composition  in  WT%  methanol 
XD  -  Overhead  Composition  in  WT%  methanol 

y  -  System  output 

y  -  Modified  system  output  defined  as  y/P 

-1  d 
z  -  Backshift  operator;  z_ku  =  u(t-k) 

t 

c.  Greek: 

e  -  Weighted  output  prediction  error 
C  -  Uncorrelated  zero-mean  random  input 

0  -  Scalar  output  function 

9  -  Vector  of  control ler  parameters:  F,  G,  H,  D 

P  -  Forgetting  factor 

o  -  Variance 

d.  Polynomials  in  z“^: 

A  -  Order  N  corresponding  to  system  output,  aQ  =  1 

B  -  Order  N  corresponding  to  control  input,  =  0 

C  -  Order  N  corresponding  to  random  disturbance, 

c  =  1 
o 

D  -  Associated  with  load  in  final  control  law 

E  -  Use  to  separate  past  and  future  disturbance 

F  -  Use  to  separate  past  and  future  disturbance 

G  -  Associated  with  control  input  in  the  final  control 
1  aw 
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H 

L 

P 

Q 

R 


Associated  with  set  point  in  the  final 

Order  N  correspondi ng  to  deterministic 
disturbance  input,  1  =  1 

Weighting  function  on  system  output 

Weighting  function  on  system  input 

Weighting  function  on  set  pont 


control 

load 


1  aw 


e.  Superscript: 


'  -  Indicates  polynomials  used  for  derivation  of  STC 
for  the  case  without  weighted  output 

~  ~  Estimates  of  the  parameters 

*  -  Predicted  value 

T  -  Transpose  of  Matrix 


f.  Subscript: 


ij  -  For  the  multivariable  case,  denotes  the  jth 
component  of  the  i th  loop 

N  -  Numerator  of  a  z  transfer  function 

D  -  Denominator  of  z  transfer  function 

t  -  Value  at  present  time 

t+k  -  Value  at  K  samples  from  present  time 


9.  Symbols: 


[  ]  -  Use  to  designate  the  variables  when  extending  the 
results  to  the  multivariable  system.  When  used 
with  polynomials  in  z  ,  the  result  is  a  poly¬ 
nomial  matrix  in  z  .  When  used  with  scalars, 
the  results  in  a  vector. 

_  -  Designates  a  vector 

_  -  Designates  a  matrix 
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APPENDIX  A 


DETAIL  SCHEMATIC  OF  PILOT  SCALE  DISTILLATION  COLUMN 


LEGEND: 


indicates  local  analog  equipment 
TYPE  -  gives  function 
TYPE.  CODES 
T  -  transmitter 
R  -  recorder 
C  -  controller 
I  -  indicator 
F  -  flow 

A  -  concentration  (analyzer) 
T  -  temperature 
P  -  pressure 
DP  -  differential  pressure 
n  -  identification  number 


HXn 


heater  exchanger 
n  -  identification  number 


gas  chromatograph 

thermocouple 

indicates  utilities 
UT  -  ST  60  psig  steam 

-  CW  -  cooling  water 


solenoid  valve 
n  -  identification  number 

control  lines 
process  lines 
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TO 


Figure  A.1  A  Detailed  Schematic  of  the  Pilot  Scale 

Distillation  Column 


APPENDIX  B 

SAMPLE  CALCULATION  TO  OBTAIN  INITIAL  Q-WEIGHTING  PARAMETERS 

1.  Using  Halman's  technique  [75],  the  kp  and  k.  are 
related  to  the  process  dynamics  according  to: 


a  = 

1  + 

0.5  *  T s / T d 

(B.4) 

k  = 

2x/3KGaTs 

(  B  .  5  ) 

ao  = 

1  + 

T$/2t 

(  B .  6  ) 

al  = 

-(1 

-  T  s  /  2  t  ) 

(  B  .  7  ) 

kp  = 

-ao 

*  k 

(  B  .  8  ) 

ki  = 

!ao 

+  3^  *  K 

(  B .  9  ) 

Where : 

kp  =  proportional  gain  constant  in  Q-wt 
k-j  =  integral  gain  constant  in  Q-wt 
Kq  =  Process  gain 

=  Process  time  delay 
Ts  =  Sample  time 
t  =  Process  time  constant 


2.  From  results  given  by  Kan  [66],  the  bottom  composition 
response  for  a  7  1/2%  decrease  in  steam  can  be  re¬ 
presented  by  a  first  order  plus  time  delay  model  with 
constants : 

Kr  =  1 . 965  wt%/ ( g/s ) 
t  =  1000  second 
T  =  180  second 
=  316  second 


3.  From  the  discussion  given  in  Chapter  2,  it  has  been 

shown  that  the  affect  of  the  process  time  delay  should 
be  removed.  So,  the  time  delay  is  set  equal  to  the 
sample  time.  Using  the  values  given  in  section  2,  the 
constants  in  section  1  become: 

a  =1.5 
k  =  1.259 
a  =1.09 
a°  =  -.91 
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The  controller  constants  are  then: 

k_  =  1.145 
kT  =  0.226 


APPENDIX  C 


PROGRAMS  LISTING 


PROGRAM  DESCRIPTION  PAGE 

STC  Implements  the  STC  Algorithm  on  the  Column  151 

HSET  Sets  Original  Controller  Parameters  and  Option  178 

ERSUM  Computes  Absolute  Error  Value  of  Data  186 

HPLOT  Plots  Data  Stored  in  Files  188 

GCLNK  Obtains  GC  Report  from  GC  Computer  196 

GCSCH  Process  GC  Report  and  Initiate  GC  Cycle  198 

GCSWT  Reset  Digital  Switch  That  Initiate  GC  Cycle  201 

BDC99  Monitors  Column  Operation  202 


0001 

0002 

0003 

0004 

0005 

0006 

0007 

0008 

0009 

0010 

001 1 

0012 

0013 

0014 

0015 

0016 

0017 

0018 

0019 

0020 

0021 

0022 

0023 

0024 

0025 

0026 

0027 

0028 

0029 

0030 

0031 

0032 

0033 

0034 

0035 

0036 

0037 

0038 

0039 

0040 

0041 

0042 

0043 

0044 

0045 

0046 

0047 

0048 

0049 

0050 
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C 

FTN4 

PROGRAM  STC  (  ,94)  ,ST  CONTROL  MOD.  BY  HYL  NOV.  1,  1979 
COMMON  I  COM ( 1 ) 

DIMENSION  IDCB2 ( 144)  ,  IDCBK144)  ,  IFILE(3)  ,  IFILU(3) 
DIMENSION  ICX(5)  ,  IDX ( 6)  ,  ERRK2)  ,  ERR2<2)  ,  IF’RAM(2) 
DIMENSION  PERERR(2),YSTAR(2),USEND(2) 

DIMENSION  YERR(2),DPR<2),BUF(56),ACLQSS(2) 

DIMENSION  NSELF(2),N(2,6),ITYID<2) , SAVE (2,1 1,30) 
DIMENSION  UN(2),DN(2),YN(2),U0(2),U0LD(2),PHI(2) 
DIMENSION  K1(2),K2(2),K3(2),UN(2),CA(2),CB(2),CQ(2) 
DIMENSION  NQ(2) ,QD(2,4) ,QN(2,4) 

DIMENSION  NP(2),PN(2,4),PD(2,4) 

DIMENSION  ITYU( 2) , IMU( 2) ,HU( 2) ,LU< 2) r BU< 2) 

DIMENSION  ITYD(2),IMD(2),HD<2),LD(2),BD(2) 

DIMENSION  ITYCT(2),UNT(2),UBL(2),UTL(2), IRAMP(2) ,ACT (2) 
DIMENSION  USS ( 2 )  ,  UPID(2) 

DIMENSION  UPROP < 2 ) ,UINT(2),UDINT(2),USINT(2),UDER(2) 
DIMENSION  AKP(2)  , AKI (2)  ,AKD(2)  ,AKU(2)  ,AKDN(2) 
DIMENSION  AKPIS(2),AKPD(2),TS(2),UDIS(2) 

DIMENSION  AKPS(2),AKIS(2),AKDS(2),ER1(2),ER2(2),SIAE(2) 
DIMENSION  ADCC(2),ADCM(2),DACC(2),DACM(2),ITEST(2,8) 
DIMENSION  KRS ( 1 6)  ,  CHAR ( 5 ) ,  ITIME ( 5) ,  I0UTCH<2) 
DIMENSION  X(2,30),PAR(2,30),P(2,210),G(2,30) ,R0U(2 ,3 ) 
DIMENSION  IPIS(2),IP1(2),IP2(2),IFP(2,2),SADC(8) 

C 

C 

DATA  NPARN/30/,NPMN/210/,NSCAN/5/,ISCAN/500/, 

*  CA  ,  CB  ,  CO  /  2*1.  ,  4*0.0/, 

*  I C  X/ 15,10,13,9,12/, 

*  ILG ,  IFILE/1 10,2H$H,2HDA,2HTA/, 

*  ITEST/1 6*0/,I0UTCH/7,8/, 

*  CHAR/4HINC  , 4HSTQ  ,4HSTI  ,4HP0S  ,4HMAN  / 

C 

C 

CALL  RMPAR(IPRAM) 

IUIQ=IPRAM( 1 ) 

I R 1 0  =  IUIO 
C 

C - - - 

c- 

-  OPEN  UP  FILE  CONTAINING  "ORIGINAL  DATA". 


C 

C- 

C- 

C- 


ICR=1 36 

5  CALL  0PENUDCB1 ,  IERR , IFILE , 0,77, ICR,  144) 


INITIALIZING  CONSTANTS. 


0051 

0052 

0053 

0054 

0055 

0056 

0057 

0058 

0059 

0060 

0061 

0062 

0063 

0064 

0065 

0066 

0067 

0068 

0069 

0070 

0071 

0072 

0073 

0074 

0075 

0076 

0077 

0078 

0079 

0080 

0081 

0082 

0083 

0084 

0085 

0086 

0087 

0088 

0089 

0090 

0091 

0092 

0093 

0094 

0095 

0096 

0097 

0098 

0099 

0100 

0101 
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C- 

C 

H=0 

DO  25  1=1,2 
C 

c 

AKDN( I )=0.0 
ERKI)  =0.0 
IPIS(I)=1 

ER2(I)  =0.0 
ERRKI)  =0.0 
ERR2 ( I )  =0.0 

S I AE ( I )  =0 . 0 
U I NT ( I  )=0.0 
US I NT ( I )=0.0 
YN( I )  =0.0 
UN ( I )  =0.0 
UNT ( I )  =0 . 0 
YS  T  AR ( I)=0.0 
UOLD(  I  )=0.0 
ACLOSS(I)=0.0 
AKNEU=0.0 
M=0 
I  JK=0 

ITEST ( I ,8)=5 
ITEST (1,1 )=1 
C 

DO  15  J=1  ,NPMN 
P(I,J)=0.0 
15  CONTINUE 

C 

DO  20  J=1 ,NPARN 
X<I,J)=0.O 
G( I , J)=0.0 
DO  20  K=1 ,11 

SAVEd  ,K,  J)=0.0 
20  CONTINUE 

25  CONTINUE 

UR  I TE ( IUI0,800) 

C 

C- - - 

c- 

C-  ASK  FOR  NO.  OF  SYSTEM  AND  INTI AL  FREQUENCY  OF  DUMPING 
C-  INFORMATION  TO  DISK. 

C- 

c 

READ  ( IRIO, *)  NSYSRD  ,  IDMP  ,  IDPR 
DO  90  IS=1, NSYSRD 
C 

C - - — - - 

c- 


0102 

0103 

0104 

0105 

0106 

0107 

0108 

0109 

0110 

0111 

0112 

0113 

0114 

0115 

0116 

0117 

0118 

0119 

0120 

0121 

0122 

0123 

0124 

0125 

0126 

0127 

0128 

0129 

0130 

0131 

0132 

0133 

0134 

0135 

0136 

0137 

0138 

0139 

0140 

0141 

0142 

0143 

0144 

0145 

0146 

0147 

0148 

0149 

0150 

0151 

0152 
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C-  READING  FROM  SECTOR  1  CIS=U  AND  4  CIS=23. 

C- 

c 

CALL  READF(IDCB1 ,  IERR,BUF,80, ILG) 

TS ( IS  )  =  BUF (2) 

CALL  READF ( IDCB1 , I ERR, BUF, 90, IDUM) 

ADCM( IS)=  BUF (40) 

ADCC( IS)  =  BUF ( 41 ) 

C 

C 

DACM( IS)  =  BUF (42 ) 

DACC( IS)=  BUF(43) 

C 

C - 

c- 

c-  READ  FROM  SECTOR  2  AND  5. 

C- 

c 

CALL  READF ( IDCB1 ,IERR, BUF, 80, ILG) 

C 

c-  NO.  OF  DELAYS  FOR  CONTROL,  INTERACTIVE  AND  FEEDFORUARD 
C-  TERM. 

C 

K 1  (IS)  =  IF I X ( BUF ( 6 ) ) 

K  2 ( IS)  =  IFIX1BUF (7) ) 

K3( IS)  =  IFIX(BUF (8) ) 

C 

C-  OBTAIN  INITIAL  PARAMETER  GUESS. 

C 

IFLG  =  0 
NSELF ( I S )  =  0 
IJ  =0 

DO  45  1=1,5 

IF <  I  .LT.  5)  N( IS, I )  =  IFIX(BUF( I) ) 

NSELF ( IS)=NSELF ( IS) +N( IS, I) 

NLOOP  =  N(ISfI) 

IF (  NLOOP  .EQ.  0)  GO  TO  45 
DO  40  J=1, NLOOP 

IJ  =  IJ  +  1 
IK  =  1*10  +  J  -  1 

IF( IK  .GT.  40)  IK  =  IK  -40 
PAR ( IS , I J )=BUF ( IK ) 

IFU.NE.4  .OR.  J.GT.1)  GOTO  40 
N ( I S , 5 >  =  I F I X ( BUF < 5) ) 

CALL  READF (IDCB1,IERR, BUF, 90, IDUM) 

IFLG  =  1 

40  CONTINUE 

45  CONTINUE 

C 

C-  TYPE  OF  IDENTIFICATION  SCHEME. 

C 


• 

0153 

0154 

0155 

0156 

0157 

0158 

0159 

0160 

0161 

0162 

0163 

0164 

0165 

0166 

0167 

0168 

0169 

0170 

0171 

0172 

0173 

0174 

0175 

0176 

0177 

0178 

0179 

0180 

0181 

0182 

0183 

0184 

0185 

0186 

0187 

0188 

0189 

0190 

0191 

0192 

0193 

0194 

0195 

0196 

0197 

0198 

0199 

0200 

0201 

0202 

0203 
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IF ( IFLG.NE. 1 )  CALL  READFUDCB1  f  IERR,BUF,90, IDUM) 
ITYID(IS)=IFIX<BUF (20) ) 

R0U< IS, 1 )=  BUF (21 ) 

R0U( IS,2)=  BUF ( 22) 

R0U( IS,3)=  BUF ( 26 ) 

PINT  =  BUF ( 23 ) 

PINTI=PINT 


IF( ITYID(IS) .NE.2)  GO  TO  50 
C 

C-  USING  SQUARE  ROOT  ID  METHOD. 

C 

c 

c 

PINT=SQRT (PINT) 

ROU( IS, 1 )=SQRT(ROU(IS, 1 ) ) 
R0U(IS,2)=SQRT(R0U(IS,2>) 
R0U(IS,3)=SQRT(R0U(IS,3)) 

50  CONTINUE 

IFP( IS, 1 )=IFIX( BUF(24) ) 

IFP< IS,2)=IFIX(BUF (25) ) 

IF ( I T Y I D ( IS)  .EQ.  3)  GO  TO  60 
C 

C-  COVARIANCE  MATRIX  FOR  ID  TYPE  1,2,  3>  4. 

C 

IJ  =  0 
IJP=  1 

NLOQP=NSELF ( IS) 

DO  55  Is  1 ,NLOOP 
IJ  =  IJ  +  I 

IF(I.NE.IFP(IS,1 ).AND.I.NE.IFP(IS,2>) 

*  P(IS,IJ)  =  PINT 

55  CONTINUE 

GO  TO  70 
C 

C-  COVARIANCE  MATRIX  INITIALIZATION  FOR  I DENT  I F I AT I ON 
C-  METHOD  3. 

C 

60  IJ  =  1 

NLOQP=NSELF ( IS ) 

IK  =  NLOOP 
DO  65  1  =  1, NLOOP 

IF(I.NE.IFP<IS,1 ).AND.I.NE.IFP(IS,2)> 

*  P(IS,IJ)  =  PINT 
IJ  =  IK  +  1 

IK  =  IJ  +  NLOOP  -  I  -  1 
65  CONTINUE 

70  CONTINUE 
C 

c- 
c 


INFORMATION  ON  SETPOINT 


. 


0204 

0205 

0206 

0207 

0208 

0209 

0210 

0211 

0212 

0213 

0214 

0215 

0216 

0217 

0218 

0219 

0220 

0221 

0222 

0223 

0224 

0225 

0226 

0227 

0228 

0229 

0230 

0231 

0232 

0233 

0234 

0235 

0236 

0237 

0238 

0239 

0240 

0241 

0242 

0243 

0244 

0245 

0246 

0247 

0248 

0249 

0250 

0251 

0252 

0253 

0254 
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c 

c- 

c 


ITYU( IS)=IFIX(BUF(30) ) 
HU ( IS)  =  BUF ( 3 1 ) 

LU ( IS)  =IFIX(BUF<32> ) 
B  U ( IS)  =  BUF(35) 

A  K  U ( IS)  =  BU< IS ) 

INFORMATION  ON  LOAD. 


C 

C- 

C- 

C- 

C 

C 

C- 

C 


C 


C 


C 

C- 

C 


ITYD(IS)=IFIX(BUF(40)) 
H  D ( IS)  =  BUF  <41 ) 

LD ( IS)  = I F I X  <  BUF ( 42 ) ) 
BD(IS)  =  BUF(45) 


READING  FROM  SECTOR  3  CIS=1]  AND  6  CIS=23. 


CALL  READFdDCBI  ,IERR,BUF,80,ILG) 
ITYCT (IS)=IFIX(BUF ( 1 ) ) 

IFdTYCT ( IS)  .NE.  3)  GO  TO  80 
DO  75  11=1,6 
N(IS,II)=0 
P  A  R  < IS, II )=0.0 
75  CONTINUE 

80  CONTINUE 


UBLdS) 
UTLdS) 
IRAMPdS) 
ACT (IS) 
NQ(IS) 

NP (IS) 


=  BUF (2 ) 

=  BUF (3) 

= I F I X ( BUF (4)) 

=  BUF(6) 

=IFIX(BUF< 1 0) ) 
=IFIX(BUF(30) ) 


DO  85  1=1,4 

QDdS,  I  )=BUF(  1+1 0 ) 

QN( IS, I )=BUF ( 1+14) 

PD( IS, I )=BUF ( 1+30) 

PN( IS, I )=BUF ( 1+34) 

85  CONTINUE 

KP,KI ,KD  USED  FOR  PID  AND  Q-UEIGHTING. 

CALL  READFdDCBI  PIERR, BUF, 90, IDUH) 
AKPS( IS)=BUF(21 ) 

AKIS(IS) =BUF(22 ) 
AKDS(IS)=BUF(23) 

AKP ( IS)=BUF (24) 

AKI(IS)=BUF(25> 


. 


. 

t 


0255 

0256 

0257 

0258 

0259 

0260 

0261 

0262 

0263 

0264 

0265 

0266 

0267 

0268 

0269 

0270 

0271 

0272 

0273 

0274 

0275 

0276 

0277 

0278 

0279 

0280 

0281 

0282 

0283 

0284 

0285 

0286 

0287 

0288 

0289 

0290 

0291 

0292 

0293 

0294 

0295 

0296 

0297 

0298 

0299 

0300 

0301 

0302 

0303 

0304 

0305 
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AKD  ( IS )  =BUF (26) 

90  CONTINUE 
C 

C-  READ  IN  PRESENT  OUTPUT  AND  EQUATE  USS  TO  PRESENT  OUTPUT. 
C 

DO  88  IS  =  1 f NSYSRD 
IDUM  =  13 

IF (IS  .EQ.  2)  IDUM  =  9 
CALL  A I RD( 5, IDUM, I  HE A, IERR) 

USS(IS)  =  10.  *  SORT ( IHEA*ADCM( IS)  +  ADCC(IS)) 

IF (ACT (IS)  .LT.  0.)  USS(IS)  =  100  -  USS(IS) 

UN ( IS)  =  USS(IS) 

UP  I D ( IS)  =  USS(IS) 

88  CONTINUE 
C 

C-  CLOSE  FILES  CONTAINING  INITIAL  STARTING  DATA. 

C 

CALL  CLOSE< IDCB1 , IERR) 

C 

c 

c 

c 

c - 

c- 

C-  ASK  FOR  NAMES  OF  TUO  FILES  UHERE  DATA  UILL  BE  DUMPED. 

C-  OPEN  FILE  ACCORDINLY. 

C- 

c 

91  U  R I T  E ( IUI0,820) 

READ(IRI0,825)IFILU,ICR 

CALL  OPEN < I DCB 1 , IERR, IF ILU, 0,77, ICR, 144) 

IF (IERR.LT.O)  CALL  FNGER( IERR) 

URITE( IUI0,820) 

READ(IRI0,825)IFILU,ICR 

CALL  0PEN(IDCB2,IERR,IFILU,0,7?,ICR,  144) 

IF (IERR.LT.O)  CALL  FMGER(IERR) 

IDISK=IDMP 
ICNT=IDISK-1 
DO  94  1=1,56 
94  BUF ( I ) =0 . 0 

MULT  =  ( (IFIX(TS(1 )) )/NSCAN) 

MULT2  =  IFIX(TS(2)/TS(1 )) 

C 

C*  IJK  =  MULT-1  UILL  CAUSE  EXECUTION  OF  CONTROL  SECTION  ON 
C*  FIRST  SCAN. 

C 

I JK=MULT-1 
I JK2  =  MULT2  -1 
C 

C-  ASK  FOR  DESCRIPTOR  FOR  EACH  DATAFILE.  ONE  LINE  OF  SO 
C-  ASCII  CHARACTERS  CAN  BE  USED. 


•] 


,*j  m  m 


0306 

0307 

0308 

030? 

0310 

0311 

0312 

0313 

0314 

0315 

0316 

0317 

0318 

0319 

0320 

0321 

0322 

0323 

0324 

0325 

0326 

0327 

0328 

0329 

0330 

0331 

0332 

0333 

0334 

0335 

0336 

0337 

0338 

0339 

0340 

0341 

0342 

0343 

0344 

0345 

0346 

0347 

0348 

0349 

0350 

0351 

0352 

0353 

0354 

0355 

0356 
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C 

DO  940  IS  =  1 , NSYSRD 
URITE(IUI0,932)  IS 
READ(IUI0,934)  (BUF(I)  ,1  =  1,20) 

IF ( IS  .EQ.  2)  GOTO  935 
CALL  URITF(IDCB1 , I ERR , BUF , 40 ) 

GOTO  937 

935  CALL  URITF < IDCB2, IERR,BUF,40) 

937  CONTINUE 

940  CONTINUE 
C 

934  FORMAT (20A4) 

932  FORMAT ( /  ENTER  DESCRIPTOR  FOR', 13) 

C 

C*  INITIALIZE  THE  START  TIME 

C*  ITM1  IS  IN  UNITS  OF  1 0*MILLISECOND. 

C 

CALL  EXEC( 1 1 ,  ITIME) 

ITM1=100*ITIME(2)+ITIME(1) 

C 

C 

c 

C*** ********************* ************************************ 

c* 

c 

c 

C*  START  OF  EACH  SCAN 
C* 

c***************************************** ************* ****** 

c 

c 

c - 

c- 

c-  READ  STATUS  OF  REGISTER  USING  SUBROUTINE  RSU.  THE  STATUS 
C-  WILL  BE  STORED  IN  ARRAY  KRS.  KRS ( 1 )  WILL  BE  THE 
C-  STATUS  OF  BIT  1  FROM  FROM  THE  LEFT, 

C-  I.E.  BIT  15  ASINDICATED  ON  THE  SUITCH. 

C-  A  VALUE  OF  1  INDICATES  THAT  THE  SUITCH  IS  UP,  0  IS  FOR  DOUN. 
C- 

c 

95  CONTINUE 
C 

CALL  EXECdl,  ITIME) 

ITM=100*ITIME(2)+ITIME(1) 

CALL  RSU(KRS) 

IF( KRS ( 1 5 ) .EQ. 1 )  UR I TE < IUIG, 1 0OO)  ITM,ITH1 , IUAIT 
1000  FORMAT*'  TIME  IS  =  ",218,'  IUAIT  =  ',14) 

C 

C-  DECIDE  UHICH  SYSTEM 
C 

IS*1 


. 


oib; 

0358 

0359 

0360 

0361 

0362 

0363 

0364 

0365 

0366 

0367 

0368 

0369 

0370 

0371 

0372 

0373 

0374 

0375 

0376 

0377 

0378 

0379 

0380 

0381 

0382 

0383 

0384 

0385 

0386 

0387 

0388 

0389 

0390 

0391 

0392 

0393 

0394 

0395 

0396 

0397 

0398 

0399 

0400 

0401 

0402 

0403 

0404 

0405 

0406 

0407 
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lHKRS(y).EU.I)  iS  =  2 
IF(KRSd)  .EQ.  1)  GO  TO  105 
IF(KRS(2)  .EQ.  0)  GO  TO  105 
C 

C - - - 

c- 

C-  BIT  2  FROM  LEFT  IS  ON,  UANT  TO  PRINT  OUT  VALUES  OF 
C-  SPECIFIED  CONSTANTS.  (BIT2  FROM  LEFT  IS  BIT  14  ON  SUITCH 
C-  REGISTER). 

C- 

C 

ASIG=+1 .0 

IF ( KRS ( 5 ) .EQ. 1 )  ASIG=-1 .0 

AKNEU=AKNEU+ASIG*(FLOAT(KRS(6))+O.1 *FLOAT ( KRS ( 7 ) )+ 

X  0.01 ♦FLOAT (KRS(8) )  ) 

IF ( KRS ( 4 ) .EQ.  1)  AKNEU=0.0 

URITE(IUI0,600)IS,AKNEU,AKP(IS),AKI(IS),AKD(IS) 

S  , AKU(IS) ,AKDN( IS) , AKPS( IS) ,AKIS(IS),AKDS(IS) 

GO  TO  110 
C 

C - - - - - 

c- 

C-  CONVERT  REGISTER  READING  TO  TYPE  OF  ITEST  AND  THE  NEU 
C-  VALUE. 

C 

c 

c- 

c 

c 

c 

105  I F ( KRS ( 1 ) .EQ.O)  GO  TO  110 

IT  =  KRS(5)+KRS(4)*2+KRS(3)*4+1 
ITEST(IS,IT)=KRS(8) +KRS(7 )*2+KRS( 6) *4 
IF (ITEST ( IS, 1 ) .EQ.7)  URITE ( IUI0,81 5 )  IDX 
C 

C - - - 

c- 

C-  CHANGE  VALUE  OF  CONSTANTS  TO  VALUE  STORED  IN  AKNEU. 

C- 

c 

110  IF(KRS(10).EQ.1)URITE(IUI0,610)((ITEST<JJ,J),J=1 ,8), 

*  JJ=1 ,2) 

IF (ITEST ( IS,4)  .EQ.  1)  AKP(IS)  =AKNEU  ♦  AKP(IS) 

IF  <  ITEST  ( IS ,  4 )  .EQ.  2)  AKKIS)  =AKNEU  +  AKI(IS) 

IF < ITEST ( IS , 4 )  .EQ.  3)  AKD < IS)  =AKNEU  +  AKD(IS) 

IF ( ITEST ( IS , 4 )  .EQ.  5)  AKPS ( IS ) =AKNEU  +  AKPS(IS) 

IF ( ITEST ( IS , 4 )  .EQ.  6)  AKIS ( IS ) =AKNEU  ♦  AKIS(IS) 

I F ( I  TEST ( I S ,  4 )  .EQ.  7)  AKDS( IS) =AKNEU  +  AKDS(IS) 

C 

IF ( ITEST ( IS,5)  .EQ.  1)  AKU( IS ) =AKHEU 
IF (ITEST ( IS,5)  .EQ.  2)  AKDN(IS)=AKNEU 
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0435 

0436 

0437 

0438 

0439 

0440 

0441 

0442 
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0444 

0445 

0446 

0447 

0448 

0449 

0450 
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0453 

0454 

0455 
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IF ( ITEST ( IS, 5)  .EQ.  3)  TPFIL=AKNEU 

IF ( ITEST ( IS , 5 )  .EQ.  6)  USS(IS)=AKNEU  ♦  USS(IS) 

IF < ITEST ( IS, 5)  .EQ.  7)  I R AMP ( I S )  =  I F I X ( AKNEU )  +  IRAHP ( I S ) 
IF ( ITEST ( IS, 6)  .EQ.  7)  MULT=( ( IFIX(AKNEU) )/NSCAN) 

C 

C-  CHECK  IF  JUST  UAITING,  UANT  TO  STOP,  READ  INITIAL  DATA 
C-  FILE  AGAIN. 

C 


IF(ITEST(1,8).EQ.5.0R.ITEST(2,8).EQ.5)  GOTO  301 
IF ( ITEST ( IS , 8 ) .GE.6 )  GOTO  999 
C 
C 

C - - - - - 

C- 

C-  CHECK  FOR  NEU  VALUE  OF  SETPOINT. 

C- 

c 


DO  157  IS=1,NSYSRD 
U0( I S )  =  BU ( IS) 

IF ( ITEST  ( IS,2) 
IF( ITEST ( IS,2) 
IF ( ITEST ( IS,2) 
IF ( ITEST ( IS, 2) 
140  IF ( ITEST ( IS ,  4 ) 

C 

C - - — 

c- 

c-  RESETTING  THE  Q-MATRIX 

C- 

C 


m 

o 

. 

0) 

I MU  < IS ) 

=0 

EQ. 

3) 

UO(IS) 

= AKU (IS) 

EQ. 

6) 

UO(IS) 

=BU( IS)+1 1 

.67 

EQ. 

7) 

U0< IS) 

=BU< IS)-1 1 

.67 

NE. 

4) 

GO  TO 

156 

NQ( IS ) =1 
DO  145  1=1,4 
C 

c 

QD( IS, I )=0 .0 
QN( IS, I )=0.0 
145  CONTINUE 


C 

c-  IF  KRS( 1 4)  IS  UP,  THEN  UANT  TO  REMOVE  Q-FILTER. 
C 

IF ( KRS ( 1 4 ) - EQ . 0 )  GO  TO  147 
QD( IS, 1 )=t 
GO  TO  155 


C 

c-  RESETTING  Q-MATRIX  TO  PID  FILTER. 

C 

147  NQ(IS)=3 

QN( IS, 1  )  =  1 .0 
QN( IS,2)=-1 .0 

QD( IS, 1 )  ?  AKP (IS)  ♦  AKI(IS)  +  AKD(IS) 
QD( IS, 2)  =  - ( AKP (IS)  +  2*AKD( IS ) ) 
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QD( IS,3)  =  AKD(IS) 

C 

155  WRITE (IUIO, 605) IS, AKP (IS) ,AKI (IS) ,  AKD( IS) 

S  ,IS,NQ(IS),(QDUS, J),J=1,4),(QN(IS,J),J=1 ,4) 

156  CONTINUE 
C 
c 

IF ( ITEST ( IS,5) .EQ.O)  GO  TO  157 

IFUTEST (IS, 5)  .EQ.4)  WRITE ( IUIO, 606)TPFIL 
IF < ITEST ( IS,5) .NE.5)  GO  TO  157 
C 

C - - - 

c- 

C-  RESETTING  THE  P-WEIGHTING. 

C- 

C 

NP( IS)  =  1 
PD( IS,  1  )  =  1 .0 
PN( IS,  1 )  =  1 .0 
PN( IS,2)=0.0 

IF(TPFIL.LE.O.O)  GO  TO  157 
NP( IS)=2 

PN(IS,2)=-EXP(-1.0/TPFIL) 

PD ( IS, 1 ) =1 .0-PN( IS, 2) 

WRITE(IWI0,607)TPFIL,IS,NP(IS),(PN(IS,J),J=1,4), 

3  (PD(IS, J) , J=1 ,4) 

157  UN( I S ) =  W  0 ( IS) 

C 

C- - - - - - - 

c- 

c-  RESETTING  THE  COVARIANCE  MATRIX. 

C- 

c 

IS  =  1 

IF ( KRS < 9 ) .EQ. 1 )  IS=2 
IF (ITEST ( IS, 6 )  .NE.3)  GOTO  480 
C 

c 

PIN=PINTI 

PIN=SQRT ( ( 1 0.0)**( AKNEW ) ) 

C 

DO  455  1=1 ,NPMN 
P(IS,I)=0.0 

IF< I -LE.30)  G(IS,I)=0.0 
455  CONTINUE 
C 

IJ=0 

NLOOP=NSELF (IS) 

DO  460  1=1 , NLOOP 
IJ=IJ+I 

IF(I.NE. IFP( IS,1 ).AND.I .NE. IFP( IS,2) )  P(IS,IJ)=PIN 
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C 

URITE(IUI0,990)IS,I,IJ,PdS,IJ) 

460  CONTINUE 
C 

480  CONTINUE 
I JK=I JK+ 1 

IF(IJK.LT.HULT)  GOTO  301 
C 

C - - 

c- 

C-  TIHE  TO  DO  CONTROL  CALCULATION. 

C-  NOTE,  HAVE  TO  CHECK  IF  DOING  BOTH  SYSTEM  OR  JUST  ONE 

C-  FOR  CASE  UHERE  2  SYSTEM  HAVE  DIFFERENT  SAMPLE  RATE. 
C-  SECOND  MUST  BE  A  MULTIPLE  OF  FIRST. 

C- 

c 


97  I JK2- I JK2+ 1 
NSYS=NSYSRD 

IF ( I JK2.LT. MULT2)  NSYS  =  1 

IF ( I JK2  .GE.  MULT  2 )  I P IS ( 2 )  =  IP  IS ( 2 >d 

IF ( I JK2.  GE.  MULT2)  IJK2  =  0 


I  JK  =  0 
M=M+1 

IF  < M . GT . 9999 )  M=1 
IPIS(I)  =  IPIS(I)  +  1 
IF (  IPISd)  .GT.  30  )  IPISd)  = 

IF (  I P I S ( 2 )  .GT.  30  )  IP  I S ( 2 )  =  1 

C 


c- 

GET  INPUT 

MEASUREMENT. 

c- 

READING  IN 

5 

POINTS,  THEY  ARE: 

c- 

POINT 

NO. 

DESCRIPTION 

c- 

15 

TOP  COMPOSITION 

c- 

10 

DISTILLATE  FLOU 

c- 

13 

REFLUX  FLOU 

c- 

9 

STEAH  FLOU 

c- 

12 

FEED 

c 

CALL  AIRD(5,ICX,IDX,IERR) 

DO  96  KJI=1 ,5 

SADC(KJI )=FLOAT ( IDX(KJI) ) 


1 


C 

c 

96  CONTINUE 

YBTM=FLOAT( I  COM < I COM <  8 )  + 1 )) +0.001 
YTOP=ADCMd  )*SADC(1)+ADCCd ) 

YNd  )  =  YTOP 
YN(2)SYBTM 
IDX(6)=YBTM 
C 

FEED  ■  ADCMd  )  *  SADC ( 5 )  ♦  ADCCd) 
DIST  *  ADCMd)  ♦  S ADC ( 2 )  ♦  ADCCd) 


r  ‘  i 


0561 

0562 

0563 

C 

r  _ 

UDISd)  =  ADCMd)  *  SADC (3 )  ♦  ADCCd ) 

U D I S ( 2 )  =  ADCMd)  *  SADC(4)  +  ADCC ( 1 ) 

U  J  0  ^ 

L 

0565 

c- 

0566 

c- 

FILTERING  OF  Y  USING  THE  P-UEIGHT . 

0567 

c- 

SUITCH  3  UP  THEN  WANT  TO  ADD  DISTURBANCE  INCREMENTALLY 

0568 

c- 

0569 

c 

0570 

DO 

159  IS=1,NSYS 

0571 

IP= 

IPIS(IS) 

0572 

D N <  IS)  *  0.0 

0573 

I F ( ITEST (IS, 3)  .EQ.  4)  DN(IS)=FEED~BD< IS)*KRS(1 

0574 

DPR(IS)=(YN(IS)-YSTAR(IS))*IDPR 

0575 

SAVE( IS, 1 ,IP)=0.0 

0576 

SAVE<IS,3,IP)=DN(IS) 

0577 

SAUE(IS,4,IP)=UN(IS) 

0578 

SAVE(IS,6,IP)=YN(IS) 

0579 

SAVE( IS,7, IP)=DPR< IS) 

0580 

159  CONTINUE 

0581 

c 

0582 

c 

0583 

DO 

240  IS=1 ,NSYS 

0584 

IP* 

IPIS(IS) 

0585 

YPD  *  YN(IS) 

0586 

NLOOP=NP( IS) 

0587 

IF ( NLOOP  .LT.  2  )  GO  TO  165 

0588 

c 

X 

0589 

c- 

YD  = 

1  /  P D ( Z )  *  YN 

0590 

c 

0591 

DO  160  J=2, NLOOP 

0592 

IC=IP- J+1 

0593 

IF ( IC  .LT.  1)  IC=IC+30 

0594 

YPD  *  YPD  -  PD<IS,J)*SA0E<IS,2,IC> 

0595 

160 

CONTINUE 

0596 

IF(PD(IS,1 ) .LE.O.O)  URITE(IUIO,700) 

0597 

IF (PD(IS,1 ) .LE.O.O)  PD< IS,  1 )  =  1 .0 

0598 

165 

YPD  *  YPD/PD < IS , t ) 

0599 

SAVE(IS,2,IP)=YPD 

0600 

c 

0601 

c- 

-  YP  * 

PN( Z)  /  PD(Z)  *  YN 

0602 

c 

0603 

YP  =  0.0 

0604 

SAVE(IS,8,IP)=0.0 

0605 

c 

0606 

c 

0607 

NLOOP=NP ( IS) 

0608 

DO  170  J*1, NLOOP 

0609 

IC=IP-J+1 

0610 

IF ( IC  .LT.  1)  IC=IC+30 

0611 

YP  =  YP  ♦  PN( IS, J)*SAVE(IS,6,IC) 

. 

( 
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YP  =  YP  -  PD(IS,J)*SAVE(IS,8,IC) 

170  CONTINUE 

YP  =  YP/PD( IS, 1 ) 

SAVE ( IS, 8,  IP)=YP 
C 

c - 

c- 

c-  ACCUMULATE  COST  FUNCTION  AND  ERROR. 

C- 

c 

IC  =  IP  -  K 1 (IS)  -  1 

IF  ( IC  .LT.  1)  IC=IC  +  30 

PHI ( IS)=  YP  -  SAVE(IS,4,IC)  +  SAVE( IS, 9, IC ) 

ACLOSS( IS)  =  ACLOSS ( IS)  +  PHI ( IS) *PHI ( IS) 

YERR( I S ) = Y N ( IS)  -  YSTAR(IS) 

C 

C  — - - - - 

c- 

C-  CHECK  UHICH  IDENTIFICATION  ALGORITHM.  AT  PRESENT  ONLY 
C-  LEAST  SQUARE  USING  LEAST  SQUARE  ROOT. 

C- 

c 

IF ( ITEST ( IS , 6 )  .EQ.  1)  GO  TO  200 

CALL  IDENT2(PHI(IS),NSELF(IS) , IS,PERERR( IS) 

*  ,X,PAR,P,G,ROU, IFP) 

200  CONTINUE 
C 


c- 

C-  CALCULATING  YSTAR  AND  X*THETA  LESS  THAT  CONTRIBUTION  OF  THE 
C-  CONTROL  ACTION  IN  THIS  SCAN. 

C- 

c 

IJ  =  0 
XP  =  0.0 
YSTAR( IS)=0.0 

c 

DO  220  1=1,5 

NLOOP  =  N( IS, I ) 

IF ( NLOOP  .EQ.  0)  GO  TO  220 
K=0 

ISTHP=IS 
ITEHP  =  I 

IF ( I  .EQ.  3)  K=K3( IS) 

IF (I.EQ.3  . AND. IDPR.EQ . 1 )  ITEMP=7 
IF ( I  .NE.  5)  GO  TO  205 
ISTMP  =  2 

IF (IS  .EQ.  2)  ISTMP  =  1 

C 

C 

ITEMP  =  1 
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K  =  K  2 ( IS) 

205  DO  215  J=1 ,NLOQP 

IJ  =  IJ  ♦  1 

IC  =  IP  -  J  +  1  -  K 

IFUC  .LT.  1)  IC=IC+30 

XP=XP+PAR(IS,IJ)*SAVE(ISTMP,ITEMP,IC) 

IF ( I . EQ - 4 )  GO  TO  210 

YSTAR( IS)=YSTAR(IS)  ♦ 

4  PAR(IS,IJ)*SAVE(ISTMP,ITEMP, IC) 

210  IC  =  IC  -  K 1 (IS) 

IF(  IC  .LT.  1)  IC  s  IC  +  30 
X(IS,IJ)=SAVE(ISTMP  f I  TEMP , IC ) 

215  CONTINUE 

220  CONTINUE 
C 

C-  CHECK  TYPE  OF  VALVE  ACTION. 

C 

IF ( ITEST ( IS, 6) .EQ.5)  ACT ( IS ) =-ACT (IS) 

IF ( ITYCT(IS)  .NE.  3)  GO  TO  225 

XP  =( YN( IS)  -  UN( IS) )  *  ABS(ACT (IS)) /ACT (IS) 
225  CONTINUE 
C 

C - - - 

c- 

C-  CALCULATING  NEU  CONTROL  ACTION. 

C- 

c- 

SAVE ( IS, 1 0, IP)=XP 
C0(IS)=0.0 

CA(IS)=QN(ISr1 )  +  PAR ( IS , 1 )+QD( IS, 1 ) 

CB( IS)=0.0 
NLOOP=NQ( IS) 

C 

DO  230  J=1 ,NLQQP 
IC=IP-J+1 

IF (  IC  .LT.  1  )  IC=IC+30 

C0(  IS)— C0(  IS)  — QD<  IS,  J)  •‘SAVE ( IS,  10f  IC) 

CO(IS)=CO(IS)-(QN(IS,J)+PAR<IS,  1  )*QD( IS, J) ) 

*  ♦SAVE ( IS , 1 , IC ) 

230  CONTINUE 
C 

c 

IF (N ( IS, 5)  .EQ.  0)  GO  TO  240 
NS=NSELF(IS)-N(IS,5)+1 
NLQQP=NQ(IS) 

CB(IS)=PAR(IS,NS)*QD(IS,1 ) 

I  J=NS 
C 

DO  235  J=1 , NLOOP 
ISTHP=2 

If(IS  .EQ.  2)  ISTMP=1 
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IF ( IC  .LT.  1)  IC=IC+30 

CQ( IS)=C0( IS)-PAR(IS,NS)*GD( IS, J) 

*  *SAVE( ISTNP, 1 f IC) 

235  CONTINUE 

240  CONTINUE 

DEN  =  CB(1)*CB<2)  -CA(1)*CA<2) 

IF (DEN  .EQ.  0.0)  WRITE ( IUIO, 705) 

IF ( DEN  .EQ.  0.0)  DEN  =1.0 
UNT (1)=(CQ(2)*CB(1 ) -CO ( 1 )+CA(2) )/DEN 
UNT(2)=(C0(1)*CB(2)-C0(2)*CA(1) )/DEN 
C 

C - - - 

c- 

C-  OBTAIN  INDEX  IP1  AND  IP2  FOR  PREVIOUS  AND  PREVIOUS, 

C-  PREVIOUS  VALUE  IN  SAVE  VECTOR. 

C- 

c 

DO  241  IS=1  ,NSYS 
IP=IPIS(IS) 

IP1 (IS)»IP-1 

I F ( I P 1 ( IS) .LT.  1 )  IP1 < I S ) =30 
I P  2  ( I S )  =  I P 1  <  I S )  - 1 
IF ( IP2( IS ) - LT . 1 )  IP2< IS )=30 
241  CONTINUE 
C 

C - - - - 

C- 

C-  CALCULATING  PID  CONTROL  ACTION. 

C- 

c 

DO  250  IS  =  1  ,NSYS 
IP=IPIS(IS) 

PCNT=(UTL( IS)-UBL( IS) )*IRAMP( IS)/1 00.0 
UTOP  =UN  < IS) +PCNT 

UBOT  =UN( IS) -PCNT 

ERR2(  IS)=  ERRKIS) 

ERR1  ( IS)=  ERKIS) 

ERMIS)  =  UN(IS)  -  YN(IS) 

C 

C-  UPDATING  SI AE  COUNTER. 

C 

IF ( ITEST ( IS,6) .EQ.6)  SIAE(IS)=O.0 

SIAE( IS )=  SIAE(IS)  ♦  ABS < ER 1 (IS) ) 
UPR0P(IS)=AKPS(IS)*ER1 (IS) 

IF ( ITEST ( IS,8) .LE. 2)  GOTO  247 

IF ( ITEST ( IS,  8) ,EQ.4.AND.AKIS(IS) .HE. 0.0)  GOTO  247 
IF ( ITEST ( IS,8) .EQ.4 )  GOTO  248 
C 
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UDER(IS)=AKDS(IS)*(ERR1 (IS)-ERI (IS) ) 
UINT(IS)=UINT(IS)+AKIS(IS)*ER1(IS) 
UPID(IS)=UPROP(IS)+UINT( IS) +UDER  < I S ) 
GOTO  246 
C 

c 

c 


247  UINT ( IS)=UN( IS)-UPROP(IS) 

IF ( ITEST ( IS, 8) .EQ.4)  GOTO  248 


C 

C - 

C- 

C-  INCREMENTAL  CONTROL  EQUATION 
C- 

c 


245  UPRQP(  IS)=  ERKIS)-ERRI(IS) 

ID1 =IP1 (IS) 

ID2=IP2< IS) 


UDER(IS)=AKDS(IS)+(YN(IS)-2.*SAVE(IS,6,ID1) 
*  +SAVE(IS,6,ID2>) 

UDINT(IS)=AKIS(IS)*ER1(IS) 

USINT ( IS)=USINT  < IS)+UDINT ( IS) 

DELTA =UPRQP(IS)+UDINT( IS)+UDER(IS) 
UPID(IS)=UPID(IS)+DELTA*AKPS(IS) 

IF ( ITEST ( IS,8)  .EQ.  4)  UPID(IS)  =  USS(IS) 

C 


246  IF ( ITEST ( IS, 8) 

IF( ITEST ( IS,8) 
IF ( ITEST ( IS  F  8 ) 

* 

IF ( ITEST ( IS,8) 
I F ( ITEST ( IS, 8) 
248  IF ( ITEST ( IS,8) 

C 

C . — 

c- 

C-  LIMIT  CHECKING. 

C- 

C 


. EQ.O )  UN(IS)=UPID(IS) 

. EQ. 3 )  UN ( IS ) =UPID( IS ) 
.NE.0.AND.ITEST(IS,8).NE.3) 

UPID(IS)  =  UN (IS) 
- EQ . 2 )  UN(IS)=UNT (IS)+USINT( IS) 
.EQ.  1 )  UN( IS)=UNT(IS) 

.EQ.4)  UN( IS) =USS (IS) 


IF (ITEST ( IS, 1 ) 
IF (UN(IS) 
IF < UP I D < IS)  .GT 
I F ( UN ( IS) 
IF ( UP I D ( IS)  .LT 


GE.  1)  U R I T E ( IUIO, 630)  UN(IS) 
GT.  UTQP  )  UN ( IS)=UTOP 

UTOP  )  UP I D < IS)=UTOP 

LT.  UBOT  )  UN( IS)=UBOT 

UBOT  )  UP  I D ( IS)=UBOT 


C 

C-  CHECK  RESET  UINDUP. 

C 

IF(UN( IS) .LT.UTL( IS) )  GOTO  249 
UINT(IS)=UINT ( IS)+(UN( IS)-UTL( IS) ) 
USINT  < IS) =US I  NT ( IS)  +  (UN( IS) -UTL ( IS) ) 
UN( IS)-UTL( IS) 
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249  iMUN(IS).BT.UBLUb))  GU1Q  Y42 

UINT(IS)=UINT(IS)+(UN( IS)-UBL( IS ) ) 

USINT(IS)=USINTdS)  +  (UN( IS) — UBL ( IS) ) 

UN(IS)sUBL(IS) 

942  CONTINUE 
C 

C . — - - - 

c- 

C-  CALCULATING  THE  NEU  YSTAR  I.E.  ESTIMATE  OF  Y 
C 

c 

c- 

ISTMP  =  2 

IF ( IS  .EQ.  2)  ISTMP  =  1 

NS  =>  NSELF(IS)  -  N(  IS, 5)  +  1 

YSTAR( IS)=YSTAR( IS) +PAR( IS, 1 )*UN( IS) 

*  +PAR<IS,NS>*UNdSTHP) 

SAVE ( IS, 1 f IP)=UN ( IS) 

C 

C-  IF  SUITCH  3  IS  UP,  UANT  TO  IDENTIFY  ONLY.  READ  IN  THE 
C-  CONTROL 
C 

IF ( K R S ( 1 4)  .EQ.  1 )  SAVE (IS, 1 f IP)  =  UDIS(IS) 

IC=IP-K1 (IS) 

IF ( IC  .LT.  1)  ICsIC+30 
X< IS, 1 )  =  SAVE  < IS, 1 ,IC) 

250  CONTINUE 
C 

c - 

c- 

C-  SENDING  OUTPUT 
C- 

c 

c 

DO  251  ISs1  ,NSYS 

USEND(IS)  =  UN ( IS) 

IF ( ACT (IS)  .LT.  0)  USEND ( IS )  =  ABS(ACTdS) )-UN(IS) 
USEND(IS)  =  USEND < IS ) **2  *  DACH(IS)  +  DACC(IS) 

IOUT  =  IF  I X  C  USEND ( IS) ) 

IF ( KRS( 1 6)  .EQ.  1)  CALL  OUT( IOUTCH( IS) , IOUT, IERR) 
IF (KRS(  1  6)  .EQ.O)  URITEdUIO,  1234) 

1234  FORMAT ( '  NOT  SENDING  OUTPUT") 

251  CONTINUE 
C 

c*** ************************************+* ******************* 

c* 

C*  SECTION  TO  PRINT  OUT  THE  RESULTS. 

C* 

c***-******* *********************** *************************** 

c 


DO  299  IS=1 , NS YS 


'  ,  M 

. 

*  y  *3  ftto 
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0867 

IP= 

IPIS(IS) 

0868 

ICH  =  I TEST  < IS,8)+1 

0869 

IF ( ITEST (IS, 1 ) .EQ.O)  60  TO  274 

0870 

IF (ITEST ( IS, 1 ) .LT.5) 

0871 

* 

URITE(IUI0,615)  CHAR( ICH) , IS, 

0872 

* 

M,UN(IS),YN(IS),YSTAR(IS),DN(IS), 

0873 

* 

UN(IS),USEND(IS),PHI(IS),SIAE(IS) 

0874 

C 

0875 

272 

IF ( ITEST ( IS, 1 ) .LT.3  .OR.  ITESTUS,  1 )  .GT.4)  GO  TO  273 

0876 

NPARsNSELF( IS) 

0877 

URITE(IUI0,625)(IS, J,PAR(IS, J) ,G(IS, J),J=1,NPAR) 

0878 

C 

0879 

273 

IF ( ITEST ( IS, 1 ) . EQ . 4 ) 

0880 

C 

0881 

C 

0882 

* 

CALL  PURITE( ITYID(IS) , 

0883 

* 

NSELF(IS),IS,IUIO,  1  ,X,PAR,P,G,RQU,IFP) 

0884 

274 

CONTINUE 

0885 

C 

0886 

r 

L 

0887 

c- 

0888 

C-  UPDATE 

THE  MEASUREMENT  VECTOR,  I.E.  ADD  EFFECT  OF  SETPOINT 

0889 

C-  AND  INTERACTIVE  TERHS. 

0890 

c- 

0891 

C 

0892 

IF (N( IS, 5) .EQ.O)  GO  TO  255 

0893 

NS=NSELF(IS)-N(IS,5)+1 

0894 

IC=IP-K2( IS)-K1 ( IS) 

0895 

IF(IC.LT.I)  I C= I C+30 

0896 

ISTMP=2 

0897 

IF ( IS.EQ.2)ISTMP=1 

0898 

X( IS,NS)=SAVE(ISTMP, 1 ,IC) 

0899 

255 

CONTINUE 

0900 

0901 

C 

p 

l - 

0902 

C- 

0903 

c-  CALCULATED  Q-FILTERED  VALUE  OF  U(T). 

0904 

c- 

0905 

C 

0906 

UQ  =  0.0 

0907 

NLOOP=NQ( IS) 

0908 

SAVE(IS,9,IP)=0.0 

0909 

IF ( ITYCT (IS)  .EQ.  3)  GO  TO  265 

0910 

DO  260  J=1 ,NLOOP 

0911 

IC=IP-J+1 

0912 

IF ( IC  .LT.  1)  IC  =  IC  ♦  30 

0913 

UQ  =  UQ  -  QD(IS,J)*SAVE(IS,9,IC) 

0914 

UQ  =  UQ  +  QN(IS,J)*SAVE(IS,1,IC) 

0915 

260 

CONTINUE 

0916 

265 

IF(QD(IS,1 )  .EQ.O.O)  UR I TE ( IUIO,710) 

0917 

IF(QD( IS, 1 ) .EQ.O.O)  QD( IS, 1 )=1 .0 

0918 

0919 

0920 

0921 

0922 

0923 

0924 

0925 

0926 

0927 

0928 

0929 

0930 

0931 

0932 

0933 

0934 

0935 

0936 

0937 

0938 

0939 

0940 

0941 

0942 

0943 

0944 

0945 

0946 

0947 

0948 

0949 

0950 

0951 

0952 

0953 

0954 

0955 

0956 

0957 

0958 

0959 

0960 

0961 

0962 

0963 

0964 

0965 

0966 

0967 

0968 
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SAVE ( IS,9fIP)=UG/QD(IS,1) 

C 

C*  ESTIMATING  U-OLD 
C 

297  IC  =  IP  -KKIS) 

I F ( I C  .LT.  1)  ICMC+30 
UOLD(IS)  *  SAVE ( IS , 4 f IC ) 

NS  =  NSELF(IS)  -  N(IS,5)  -  N(IS,4>  +  1 

NF  =  NSELF (IS)  -  N(IS,5) 

C  IF ( ITYCT (IS)  .EQ.  1)  GO  TO  290 

IF ( ITYCT (IS)  .EQ.  3)  GO  TO  290 
DO  298  I=NS , NF 

ICMP  -  I  -  KKIS)  *  NS 

IF ( IC  .LT.  1)  ICMC+30 

X( IS, I )=X( IS, I )  -  SA0E(IS,9,IC) 

298  CONTINUE 
290  CONTINUE 

C 

C 

299  CONTINUE 
C 

c 

C - - - - 

c- 

C-  CHECK  IF  IT  IS  TIME  TO  STORE  INFORMATION  TO  DISK. 
C- 

ICNTMCNT  +  1 

IFUCNT.LT. IDISK)  GOTO'  301 
C 

C - - - - - 

c- 

c-  TRANSFER  INFORMATION  TO  DISK. 

C- 

c 

I C  NT  =  0 
IDISKMDMP 
C 

DO  305  ISM  , NSYS 

IF ( ITEST ( IS,7) . E Q . 4 )  GOTO  999 

IF  ( ITEST ( IS,7)  .EQ.7 )  IDMPMFIX  ( AKU ( IS) ) 

IF ( ITEST < IS, 7) .LE .  1 )  GOTO  307 
BUF ( 1 )=FLQAT ( ICH) 

BUF(2)=FL0AT (IS) 

BUF  <  3 ) =FLO AT (M) 

BUF (4)=UN(IS) 

BUF (5)=YN( IS) 

BUF(6)=UN( IS) 

BUF(7)=PHI (IS) 

BUF (8)=SIAE( IS) 

BUF(9)=YT0P 
JUF ( 1 0)=UDIS( 1 ) 
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B  UF ( 1 1 )=YBTM 
BUF ( 1 2)=FEED 
BUF ( 13)=UDIS(2) 

BUF ( 1 4)=DIST 
BUF ( 15)=YSTAR( 1 ) 

BUF ( 1 6)=YSTAR(2) 

C 

c - - - 

c- 

C-  BUF (15)  -  BUF (85) 

C- 

c 

c 

NPAR=NSELF (IS) 

NPTS  =  2  *  <16  +  2*NPAR) 

DO  302  I=1,NPAR 

302  BUF (16+1 )=PAR( IS, I ) 

C 

DO  303  I=1,NPAR 

303  BUF ( 1 6+NPAR+I )=G( IS , I ) 

C 

C 

c 

c 

C  IF (HPAR.GE. 1 0)  GOTO  308 
C  NPE=(NPAR**2+NPAR)/2 

C  DO  304  1=1, NPE 

C304  BUF ( 1 4+2*NPAR+I )=P( IS,1) 

C 

308  IF ( IS.EQ.2)  GOTO  306 

CALL  URITF(IDCB1 ,IERR,BUF,NPTS) 

GOTO  307 

306  CALL  UR  I TF  < IDCB2, IERR,BUF,NPTS) 

307  CONTINUE 
C 

305  CONTINUE 
C 

C*  CALL  TIME  TO  CHECK  HOU  LONG  PROGRAM  SHOULD  BE  SUSPENDED 
C*  BEFORE  NEXT  SCAN.  NOTE,  EXCEPT  FOR  INITIAL  SCAN,  THE 
C*  START  TIME  OF  THE  NEXT  SCAN  IS  CALCULATED  FROM  END  OF 
C»  THIS  SCAN  +  TIME  IN  UAIT.  IF  THE  PROGRAM  BECOMES 
C*  SUSPENDED  FOR  LONGER  THAN  SPECIFIED  TIME  (DUE  TO  LOAD- 
C*  ING)  THE  NEXT  SUSPEND  UILL  BE  SHORTENED. 

C 

301  CONTINUE 

CALL  EXEC ( 1 1 , I T I ME ) 

ITM2=100*ITIME(2)+ITIME(1) 

IDIF=ITM2-ITM1 

312  IF(IDIF.GE.O)  GO  TO  315 
IDIF=IDIF*6000 
GO  TO  312 
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C 

C*  CHECK  IF  CONTROL  CALCULATION  UITHIN  THE  SCAN  TIHE. 

C 

315  IF (IDIF.LE.ISCAN)  GO  TO  320 
URITE(IUI0,9Q0)  IDIF 
IDIF=IDIF-ISCAN 
ITM1 =ITM1 +ISCAN 
I JK=I JK+1 

IF (IJK.GE.MULT)  GO  TO  97 
GO  TO  315 
C 

320  IUAIT=ISCAN-IDIF 
ITM1 =ITM2+IUAIT 
CALL  UAIT<IUAIT,0,IERR) 

GO  TO  ?5 
C 

C - 

c- 

C-  CLOSE  FILES  AND  EXIT 
C- 

c 

9??  CALL  CL0SE(IDCB1 f IERR) 

CALL  CLOSE( IDCB2, IERR) 

C 

IF ( ITEST ( IS, 8)  .EQ.  6)  GO  TO  5 
C 
C 

IF ( ITEST ( IS,7)  .EQ.  4)  GO  TO  91 
STOP 
C 

C************************************* ********* ************** 

c* 

C*  FORMAT  STATEMENTS 
C* 

c *********************************************** ************* 

c 

600  FORMAT ( "  AK-",I1 , "=",G12.5,"  P",F8.2,"  I",F8.2, 

S"  D",F8.2,"  U",F8.2  ,"  D",F8.2,/, 

S  "  PID  =",F8.2,"  I",F8.2,"  D",F8.2) 

605  FORMAT ( "  RESET  Q-FILTER  IS=",I1,/, 

S  IX/  AKP=",F8.3,"  AKI=",F8.3,"  AKD=",F8.3,/, 

S  "  QD(",I1,",',I1,")=",G12.5,1X,G12.5,1X,G12.5,1X, 

SGI  2. 5,//  QN(-,-)=",G12.5,1X,G12.5,1X,G12.5,1X,G12.5,  IX) 

606  FORMAT ( "  T-P-FILTER  =",612.5) 

607  FORMAT ( "  RESET  P-FILTER  T-P=" ,G12.5, /, 

S  "  PD(^,I1fx^,I1fx)  =  x,G12-5,1X,G12.5,1XfG12-5,1X, 

SG12.5,/,"PN(-,-)=",G12.5,1X,G12.5,1X,G12.5,1X,G12.5) 

610  FQRMAT("  ITEST=" ,8< 12, 1 X) , 5X ,8( I2f 1 X )  ) 

615  FORMAT ( 1 OX, "ITR",2X, "SETPT" ,2X,"MEAS.",3X,"YSTAR" ,5X, 
*"DN" ,6X, "UN" ,5X, "OUT" ,6X, , "PHI" ,5X,"SIAE" ,/ , 

*  IX, A4, "M-", II, "  =  ", 15, IX, F5. 2, IX, F7. 3, IX, F7. 3, IX, 
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&F/.2,1X,F/.3,1X,F7.1,1X,G8.3,1X,G8.3) 

625  FORMAT < IX, "PAR(",  II , ,  13," >  =  "  ,G1 2.5,"  K-G=",G12.5) 

630  FORMAT ( 1 X, 'CAL.  OUTPUT  *  ",G10.5) 

700  FORMAT ( "  PD  IS  ZERO,  RESET  TO  1.0  ") 

705  FORMAT ( x  CA  IS  ZERO,  RESET  TO  1.0  "> 

710  FORMAT ( "  QD  IS  ZERO,  RESET  TO  1.0  ") 

800  FORMAT (2X, "ENTER  NSYS, IDMP, IDPR:  V  ') 

805  FORMAT (217) 

810  FORMAT (IX, 5(612.5, IX)) 

815  FORMAT ( 1 X ,  "TP-PR  ",F7.2,"  DS-FL  ",F7.2,"  RF-FL  ",F7.2, 
i  "ST-FL  ',  F7.2,/  FE-FL  ",F7.2,"  BT-PR  ",F7.2) 

820  FORMAT (2X, 'INPUT  DATA  DUMP  FILES . MUST  EXIST", 

*  "  ALREADY",/, 

&  2X,"ITEST(7)  MUST  BE  2  TO  DUMP  !!!!") 

825  FORMAT (3A2,I3) 

900  F0RMAT(2X, "  EXCEED  TIME  SCAN,  TIME  NEEDED  =  ",I10> 

990  FORMAT ( "  IS=",I2,'  L00P=",I7,"  P< " , 14, " )=" ,G12.5) 

END 
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FTN4 

SUBROUTINE  ZTRAN(ITYP,TS, GAIN,T,A,B,NA,NB> 

DIMENSION  T(3),A(20),B<20) 

DO  10  1=1,20 
A ( I ) =0 . 0 
B( I )=0.0 
10  CONTINUE 
N  A = 1 
NB=1 

A(1)=1.0 
B( 1 )=0.0 

IF( ITYP  .LE.  0  .OR.  ITYP  .GT.  5)  GO  TO  40 
GO  T0(  15  ,  20  ,  25  ,  30  ,  35)  ,  ITYP 

CC  FIRST  ORDER  LAG  1/(TS+1)  . 

C 

15  NB=1 
NA=2 

B(1)=GAIN*<1.0-EXP(-TS/T<1)  )  ) 

A(2)=-EXP(-TS/T ( 1 ) ) 

GO  TO  40 
C 

C  SECOND  ORDER  1/(S)(S+1) . 

C 

20  NB=2 

NA=3 

B ( 1 )=GAIN*(TS-T(1 )*(1 . -EXP ( — T S/ T ( 1 ) )  )  ) 

B(2)  =  GAIN*<  -EXP ( -T S/T ( 1 ) ) *TS  +  T( 1 >*< 1 .-EXP(-TS/T( 1 ) )  ) 
A(2)=-(1.0+EXP(-TS/T(1))> 

A(3)=-EXP(-TS/T (1 ) ) 

GO  TO  40 
C 

C  SECOND  ORDER  1/(ST1 +1 ) (ST2+1 ) 

C 

25  AA=EXP(-TS/T ( 1 ) ) 

BB=EXP(-TS/T (2) ) 

NB=2 

NA=3 

B( 1 )=GAIN*(  (BB/T (2 )-AA/T ( 1 ) )/( 1 . /T( 2 )-1 ,/T< 1 ) ) - ( AA+BB) +1 
B(2)=GAIN*(AA*BB-(BB/T(2)-AA/T(1 ))/(1 ./T(2)-1 ,/T( 1 ) )  ) 
A(2)=-(BB+AA) 

A(3)=AA*BB 
GO  TO  40 
C 

C  SECOND  ORDER  1 / ( TS  ♦  1 ) < TS  +  1) 

C 

30  N  B  =  2 

NA=3 

AA=EXP(-TS/T ( 1 ) ) 

B(1 )=GAIN*(  -(1 . +TS/T( 1 ))*AA+1 .0) 

B ( 2)  =  GAIN*. ( AA*AA  -  (1.  -  TS/T <  1 ) ) *AA ) 

A(2)=-2.*AA 
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0052 

A(3)=AA*AA 

0053 

GO  TO  40 

0054 

C 

0055 

C 

SECOND  ORDER  1/<T1S*S  +  T2S  +  1) . 

0056 

C 

0057 

35 

AA  3  T (2)/(2.*T ( 1 ) ) 

0058 

BB  3  1 / T ( 1 >  -  (T(2)*T(2))/(T(1)*2.)**2 

0059 

IF(BB.LT.O.O)  GO  TO  40 

0060 

BB=SQRT(BB) 

0061 

B( 1 )=EXP(-AA*TS)*(-CQS(BB*TS)-AA/BB*SIN 

0062 

B(2)=EXP(-2.*AA*TS)-EXP(-AA*TS)*(C0S(BB 

0063 

S 

i  /BB*SIN( BB*TS) ) 

0064 

A( 1 )=1 .0 

0065 

A(2)3-2.*EXP(-AA*TS)*C0S(BB*TS) 

0066 

A(3)=EXP(-2.*AA*TS) 

0067 

A(3)=EXP(-2.*AA*TS) 

0068 

B(1 )=GAIN*B(1 )/(T(1) *<BB*BB+AA*AA>> 

0069 

B(2)=GAIN*B(2)/(T ( 1 )*(BB*BB+AA*AA) ) 

0070 

NB=2 

0071 

NA=3 

0072 

40 

RETURN 

0073 

END 

0179 

SUBROUTINE  RSU(IRSU) 

0180 

DIMENSION  I RSU ( 16) 

0181 

CALL  DI(1,1,ID,IERR) 

0182 

IRSU( 1 6)=0 

0183 

IF < ID  .LT.  0)  I RSU ( 1 6 )S1 

0184 

IF ( ID  .LT.  0)  ID  =  ID  +  32766 

0185 

I F ( I ERR  .NE.  1)  GO  TO  999 

0186 

DO  100  I  =  2,16 

0187 

IJ  =  16  -  I 

0188 

I RSU ( I J+1 )=0 

0189 

IDT  =  ID  /  (2**1 J) 

0190 

I F < IDT  .EQ.  1)  IRSU( IJ  +  1 )S1 

0191 

IF (IDT  .EQ.  1)  ID  3  ID  -  2* 

0192 

100 

CONTINUE 

0193 

DO  200  I  3  1,8 

0194 

ITEMP  3  IRSU (17-1) 

0195 

I RSU ( 1 7- 1 ) 3  IRSU ( I ) 

0196 

I RSU ( I )  3  ITEMP 

0197 

200 

CONTINUE 

0198 

999 

RETURN 

0199 

END 

2 
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0046 
0047 
0048  C 
0049  C 
0050  C 
0051  C 
0052  C 
0053  C 
0054  C 
0055  C 

SUBROUTINE  IDENT2(YNEU,N, IS, PERER, X, A,P,G,ROU, IFP) 
DIMENSION  X(2,30),A(2,30),P(2,465),G(2,30) ,RQU ( 2 ,3 ) , IFP (2, 

RECURSIVE  SQUARE  ROOT  ALGORITH . 

0056  C 

0057 

PERER  =  YNEU 

0058 

DO  10  1  =  1, N 

0059 

PERER  =  PERER  -  X( IS, I )*A( IS, I ) 

0060 

10 

CONTINUE 

0061  C 

0062 

GAMA  =  ROU( IS,  1 ) 

0063 

GAHA2=  ROU ( IS ,  1 )  *ROU( IS, 1 ) 

0064 

I  J=0 

0065 

J I  =0 

0066 

DO  30  J=1,N 

0067 

PX  =  0.0 

0068 

J1  =  J  -  1 

0069 

DO  15  1=1, J 

0070 

JI  =  JI  +  1 

0071 

PX  =  PX  ♦  P(IS,JI)*X(IS,I) 

0072 

15 

CONTINUE 

0073 

ALFA  =  GAHA/ROU( IS,2) 

0074 

BETA  =  PX/GAMA2 

0075 

GAMA2  =  GAMA2  '+  PX*PX 

0076 

GAMA  =  SQRT(GAMA2) 

0077 

ALFA  =  ALFA/GAMA 

0078 

G( IS, J)  =  P(IS,JI)*PX 

0079 

P  < I S , JI )  =  ALFA*P(IS, JI) 

0080 

IF ( JI  .EQ.  0)  GO  TO  25 

0081 

DO  20  1=1, JI 

0082 

IJ  =  IJ  +  1 

0083 

PQP  =  P(IS,IJ) 

0084 

P(IS,IJ)=ALFA*(PQP-BETA*6(IS,I)  ) 

0085 

G( IS,I)=G(1S,I)  +  PQP*PX 

0086 

20 

CONTINUE 

0087 

25 

CONTINUE 

0088 

IJ  =  IJ  ♦  1 

0089 

30 

CONTINUE 

0090  ( 

% 

0091 

PERER  =  PERER 

0092 

DO  35  1  =  1  ,N 

0093 

G(IS,I)=G(IS,I)/GAMA2 

0094 

A( IS,  I )  =  A( IS, I)  ♦  G(IS,I)*PERER 

0095 

35 

CONTINUE 

0096  1 

C 

0097 

RETURN 

0098 

END 

■ 


• 
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SUBROUTINE  PURITE( ITYPE.N, IS , IUIO, IPRNT,X, A,P,G,RQU, IFP) 
DIMENSION  D( 30)  ,  CHAR(5) 

DIMENSION  X<2,30),A(2,30),P(2,465),G(2,30) , R0U(2,3) , IFP (2, 2) 
C . 

c 

C  PRINTING  ROUTINE  . 

C 

C 

C 

c . 

c 

c 

DATA  CHAR/4HRLS  ,4HRSR  ,4HUDU  , 4HRL  ,4H+— / 

JS  =  0 

IFIITYPE.LT.I  .OR.  ITYPE  .GT.  3)  GO  TO  65 
IF  ( IPRNT.EQ.  1 )  URITEdUIO,  70)  CHAR  (ITYPE) 

DO  60  1*1, N 
DO  10  L=1  ,N 
10  D(L)*0.0 

GO  TO  (  15  ,  25  ,  40  ,  65), ITYPE 
C 

C . . -  RLS - 

C 

15  DO  20  J=I ,N 

IJ=  JMJ-D/2  +  I 
I F < J  .IT.  I)  IJ=  I*< 1-1 )/2  +  J 
D<  J-I  +  1 )=P( IS, I J) 

20  CONTINUE 

GO  TO  55 
C 

C- . -RSR - - 

C 

25  DO  35  J= I , N 

KDO  =  N  -  J  +  1 
J1  =  J 

DO  30  K=1 ,KDO 

J2  *  J1 *( J1-1 )/2  +  I 
J3  =  J1*( J1-1 )/2  +  J 
D( J-I+1 )=D< J-I+1 )  +  P( IS, J2)*P( IS, J3) 

J1*J1+1 


30  CONTINUE 

35  CONTINUE 

GO  TO  55 
C 

c— . —  U  -  D  -  U  — 

c 

J3  =  JS 
JDO=  N  -  I  +  1 
DO  50  J*1 , JDO 
D( J)*0.0 


40 


■ 


014? 

0150 

0151 

0152 

0153 

0154 

0155 

0156 

0157 

0158 

0159 

0160 

0161 

0162 

0163 

0164 

0165 

0166 

0167 

0168 

0169 

0170 

0171 

0172 

0173 

0174 

0175 

0176 

0177 

0178 
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SUM  =1.0 

KDO  =  JDQ  -  J  +  1 
DO  45  K=1 ,KDQ 

J1  =  JS  +  (J-1)  ♦  K 

J3  =  J3  +  1 

J2  =  J2  ♦  KDO  -  K  +  2 

I F < K  .EQ.  1)  J2  =  J3 

IF(J.GT.I)  SUM  =  SUM  *  P(IS,J1> 

IF(K.GT.I)  SUM  =  SUM  *  P(IS,J3) 

SUM  =  SUM  *  P(IS,J2) 

D(J)=  D( J)+SUM 
45  CONTINUE 

50  CONTINUE 
JS  =  J1 


-  PRINT  SECTION  - 

55  ND=  N  -  I  +  1 

URITEUUIO,  80)ND,(D(L),L=1,ND) 

IF ( ND.GT.5)ND=5 

IF ( IPRNT .EQ.  1 )  URITEUUIO,  75) <CHAR<5) ,L=1 ,ND) 

60  CONTINUE 
65  RETURN 

70  FORHAT  ( '  +-SA4/-+-  DIAGNAL 

l  /  + - ♦  - - -+') 

75  FORMAT  ( '  + . ',5(A4,9H -  )) 

80  FORMAT ( *  I  P',13,'  I',20(G12.5,'I',G12.5,'I',G12.5,'I', 
l  G12.5,'I',G12.5,'I',/,'  I  CONT  I')  ) 

END 
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FTN4 

PROGRAM  HSET  ()  ,  THIS  PROGRAMME  GENERATES  SELF-TUNER  DATAFI 
DIMENSION  D( 6,85)  ,  N(10)  ,  AIN(IO)  ,RTYPE(6) 

DIMENSION  I  DC  B  < 1 44)  ,  BUF( 45)  ,  IFILE(3) 

DIMENSION  PRIDNT (5) ,PRITYP( 3) 

DIMENSION  T(3)  ,  A(20)  ,  B ( 20 ) 

DATA  PRIDNT/4HN0NE,4HRLS  ,4HRSR  ,4HU-D  ,4HRL  / 

DATA  RTYPE/4HT0P  ,4HT0P  , 4HT0P  , 4HB0T  ,4HB0T  ,4HB0T  / 

DATA  PRITYP/4HST  , 4HSTPI , 4HPID  / 

DATA  KC,KM,KS,KR/1 HC,1 HM, 1HS, 1  HR/ 

DATA  KTU0/1H2/ 

DATA  ILP/2HLP/ 

C 

C 

DATA  ILEN, IFILE/1 70,2HXH,2HDA,2HTA/ 

C 

C  ...<  I/O  ROUTINE  > . 

C 

IU=LGGLU< IDUMMY) 

I R  =  I U 

CALL  OPEN< IDCB,IERR ,IFILE, 0,77,136, 144) 

IFUERR  .GE.  0)  GO  TO  5 
CALL  FMGER(IERR) 

STOP  4 
C 
C 

5  UR  I TE ( I U , 3  20 ) 

READUU,*)  IFLG 

320  FORMAT ( '  ENTER  1  IF  INITIALIZING  FILE:  V 
IF(IFLG.EQ.I)  GO  TO  15 
DO  10  J=1 ,6 

CALL  READF ( IDCB, IERR,BUF ,80, ILEN) 

IF  < I  ERR  .GE.  0)  GO  TO  6 
CALL  FMGER ( I  ERR ) 

STOP  4 

6  DO  800  JJ  =  1,40 

D(J,JJ)  =  BUF ( JJ) 

800  CONTINUE 
C 

CALL  READF ( IDCB, IERR, BUF, 90, ILEN) 

DO  802  JJ  5  41,85 
IDUM  =  JJ  -  40 
D(J,JJ)  =  BUF ( IDUM) 

802  CONTINUE 
10  CONTINUE 
15  URITE< IU, 1 25) 

IG=0 

READ  ( IR,310).ICHAR,ICHARS,IL0G 

IF ( ICHARS.EQ.KTUO)  IG=1 

IF ( ICHAR  .EQ.  KR  )  GO  TO  90 

IF( ICHAR.EQ.KC  .OR. ICHAR. EQ. KM  .OR. ICHAR. EQ.KS  ) 


. 
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♦  GO  TO  20 
GO  TO  15 

20  ISYS=1  +  IG*3 

IF ( ICHAR  .EQ.  KS  )  GO  TO  35 

IFdCHAR  .EQ.  KC  )  GO  TO  45 

C 

C . MODEL  SECTION . 

C 

IF ( ILOG.EQ .  ILP )  10=6 
NTYPE=IFIX(D(ISYS,1>) 

URITE( 10, 1 500)  RTYPE(ISYS) 

OR I TE ( 10,1 50 ) NT  YPE 
IF(NTYPE  .GT.  0)  GO  TO  25 
C 

C .  Z  -  TRANSFER  MODEL  SIMULATION . 

C 

NF  =IFIX(D( ISYS, 10)) 

NFP=NF+1 1 

ORITEdO, 155)  NF, (DdSYS, J),J=  11, NFP) 
ORITE< 10,120) 

NF=IFIX(D( ISYS, 20) ) 

NFP=NF+21 

URITE< 10,1 60)NF, <D< I SYS , J ) , J=2 1 ,NFP) 
OR I TE ( I U , 1 20) 

NF=IFIX(D( ISYS, 30) ) 

NFP=NF+31 

ORITE ( 10, 1 65)NF, (D  < ISYS , J ) , J=31 ,NFP) 
URITE ( 10,120) 

NF=IFIX(D(ISYS,40) ) 

NFP=NF+41 

ORITE (10,1 70 )NF,(D( ISYS, J ), J=41 ,NFP) 
ORITE( 10, 120) 

NF=IFIX(D( ISYS, 50) ) 

NFP=NF+51 

ORITE ( 10, 1 75) NF,(D( ISYS, J),J=51 ,NFP) 
ORITEdO,  120) 

GO  TO  30 
C 

C . RK  -  SIMULATION . 

C 

25  ORITEdO,  1  80)  (DdSYS,  J) ,  J  =  1 0,1  4) 
ORITEdO, 120) 

ORITE (10,1 85) (D (ISYS, J),J=20,24) 
ORITEdO, 120) 

ORITEdO, 1 90  )<D(ISYS,J) ,  J=30,34 ) 
ORITEdO, 120) 

ORITEdO,  1  95  )(D  (ISYS,  J),J=40, 44) 
ORITEdO, 120) 

30  ORITE(  10,200)  ( DdSYS,  J),J=2, 6) 

DELAY  =  DdSYS, 9)-D(  ISYS, 7) 

ORITE ( 10,205) (D( ISYS, J) , J=7,9) ,DELAY 


. 
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0103  URITEdU, 210)(DdSYS,J),J  =  80, 83) 

0104  UR  I TE ( IU, 1 20) 

0105  I U = I R 

0106  GO  TO  70 

0107  C 

0108  C . SYSTEM  SECTION . 

0109  C 

0110  35  ISYS=2  +  I G  *  3 

0111  DO  40  J=1,9 

0112  N( J)=IFIX(D( ISYS, J) ) 

0113  40  CONTINUE 

0114  IF ( ILOG.EQ. ILP)  IU=6 

0115  NFIX=N(1)+N(2)+N(3)+1 

0116  U R I TE < IU,21 50)  RTYPE(ISYS) 

0117  URITEdU,215)(N(J),J=1,8),NFIX,N<9) 

0118  URITE(IU,120) 

0119  NF=N( 1 )+10 

0120  URITE dU, 220) (D( ISYS, J),J=10,NF) 

0121  NF=N(2)+20 

0122  URITEdU, 225)  (D( ISYS,  J),J=20,NF> 

0123  NF=N(3)+30 

0124  URITE  dU,230)<  DdSYS,  J),J=30,NF) 

0125  NF  =  N(  4 ) +40 

0126  URITEdU, 235)(DdSYS,J),J= 40, NF) 

0127  NF  =  N  <  5 ) +50 

0128  URITEdU, 240)  (DdSYS,J),J=50,NF) 

0129  URITEdU, 120) 

0130  NF=IFIX(D( ISYS, 60) ) 

0131  NFP  =  NF+1 

0132  IF ( NFP  .LT.  1  -OR.  NFP  .GT.  5)  NFP  =  1 

0133  N1 =IFIX ( D(ISYS,64) ) 

0134  N2- 1 F I X  <  D ( ISYS, 65) ) 

0135  STD=2 . 75*SGRT (  ABS (  (D(ISYS-1 ,5)**2+0. 00001 J/30.0  )) 

0136  URITE (IU, 245)NF, PRIDNT ( NFP), (D (ISYS, J),J=6I ,63), N1,N2,D(ISY 

0137  $  ,D(ISYS,67) , STD 

0138  URITEdU, 120) 

0139  URITEdU, 250)  (DdSYS,  J) ,  J=70, 75) 

0140  URITEdU, 120) 

0141  URITEdU, 255)  (DdSYS,  J)  ,J  =  80, 85) 

0142  URITEdU, 120) 

0143  IU=IR 

0144  GO  TO  70 

0145  C 

0146  C . CONTROLER  SECTION . 

0147  C 

0148  45  ISYS=3  +  IG*3 

0149  N1=IFIX(D< ISYS, 1 ) ) 

0150  IF ( N1 . LT .  1  .OR.  N1.GT.3)  N1  =  1 

0151  N2=IFIX(D( ISYS, 5) ) 

0152  NQ=IFIX ( D( ISYS, 10)) 

0153  NP=IFIX(D( ISYS, 20)) 


. 

m 
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0154  N3=IFIX(D( ISYS, 30) ) 

0155  NR=IFIX(D( ISYS, 40) ) 

0156  PCNT=(D(ISYS,3)-D(ISYS,2)>*D<ISYS,4>/100.0 

0157  IF ( ILOG.EQ. ILP)  IU=6 

0158  UR I T E ( IU,2600)  RTYPE(ISYS) 

0159  UR  I TE  < IU,260)N1  ,PRITYP< N1 ) , ( D< ISYS , J ) , J=2 , 4) ,PCNT 

0160  C 

0161  C . Q-FILTER . 

0162  C 

0163  URITE(IU,265)N2,D(ISYS,6) 

0164  IF ( N2  .EQ.  0)  GO  TO  55 

0165  URITE(IU,27QHD(ISYS,J),J=7,9) 

0166  URITE( IU, 1 20) 

0167  C 

0168  C . ESTIMATION  OF  EQUILANT  Z-TRANSFORM . 

0169  C 

0170  NQ  =  1 

0171  D( ISYS, 1 5)=1 .0 

0172  D( ISYS, 1 6)=0.Q 

0173  IF ( D ( ISYS f 7 )  .GT.  0.0)  GO  TO  50 

0174  I F ( D < ISYS, 8)  .EQ.  0.0)  GO  TO  55 

0175  C 

0176  C . INTEGRAL  ACTION  ONLY  . 

0177  C 

0178  NQ=2 

0179  D( ISYS, 1 1  )  =  D( ISYS, 8) 

0180  D( ISYS, 1 2)=0.0 

0181  D( ISYS, 1 6)=-1 .0 

0182  GO  TO  55 

0183  50  NQ=1 

0184  D< ISYS, 11 )=D(ISYS,7)+D(ISYS,8) 

0185  D( ISYS, 1 2 )  =0 . 0 

0186  D( ISYS, 1 3)=D( ISYS, 9) 

0187  D( ISYS, 1 6)=0.0 

0188  I F ( D ( ISYS, 8)  .LE.  0.0)  GO  TO  55 

0189  NQ=2 

0190  D(ISYSf 12)=-D(ISYS,7)-2.*D(ISYS,9) 

0191  D( ISYS, 1 6)=-1 .0 

0192  IF ( D( ISYS, 9)  .GT.  0.0)  NQ=3 

0193  55  URITE(IU,275)NQ,(D(ISYS,J),J=11,14) 

0194  U  R I T  E ( I U , 2  80 )  < D( ISYS , J) , J=1 5, 1 8) 

0195  D ( ISYS , 1 0 ) =NQ 

0196  C 

0197  C .  P  -  FILTER  . 

0198  C 

0199  UR  I TE ( IU, 1 20) 

0200  NA  =  1 

0201  ITYP  =  I F I X ( D ( ISYS, 20) ) 

0202  TS  =  D( ISYS-2,2) 

0203  URITE( IU, 285) ITYP, TS 

0204  IF ( ITYP.NE.  1)  GO  TO  60 


•  . 


0205 

ITYPS=IFIX<D(ISYS,21)) 

0206 

GAIN  =  D(ISYS,22) 

0207 

T ( 1 )  =  D  < I  SYS ,23  > 

0208 

T ( 2 )  =  D ( ISYS,24 ) 

0209 

T(3)  =  B(ISYS,25) 

0210 

URITE(IU,290)ITYPS,GAIN,(T(J),J=1,3) 

0211 

CALL  ZTRAN(ITYPS,TS,GAIN,T,A,B,NA,NB) 

0212 

D(ISYS,30)=FLQAT(NA) 

0213 

D( ISYS,31 )=B( 1 ) 

0214 

D(ISYS,32)=B<2) 

0215 

D( ISYS,33)=B<3) 

0216 

D( ISYS,34)=B(4) 

0217 

D( ISYS,35)=A( 1 ) 

0218 

D( ISYS,36)=A(2) 

0219 

D( ISYS,37)=A<3) 

0220 

D(ISYS,38)=A(4) 

0221 

60 

URITE(IU,295)NA,<D<ISYS,J),J=31,34),(D(ISYS 

0222 

C 

0223 

C. 

R  -  FILTER . 

0224 

C 

0225 

URITE( IU, 1 20) 

0226 

URITE< IU, 300) ( D ( ISYS , J) , J=50  f  54 ) 

0227 

URITE(IU,301)D<ISYS,60) 

0228 

URITE(IU,302)(D(ISYS,J),J=61,63) 

0229 

URITE(IU,303)(D(ISYS,J),J=64,66) 

0230 

65 

I SY  S  =  I  SYS 

0231 

URITE( IU, 1 20) 

0232 

I U = I R 

0233 

C. 

..INPUT  OUTPUT  SECTION 

0234 

c 

0235 

70 

URITEdU, 130)1  CHAR, I  CHARS,IG 

0236 

READ  ( IR,31 0) ICHAR2 

0237 

c 

0238 

c. 

..CHANGE  ,  CHANGES  ,  RESET  ,  PRINT  . . 

0239 

c 

0240 

IF ( ICHAR2  .EQ.  KC  )  GO  TO  75 

0241 

IF ( ICHAR2  .EQ.  KS  )  GO  TO  80 

0242 

IF ( ICHAR2  .EQ.  KR  )  GO  TO  15 

0243 

GO  TO  70 

0244 

c 

0245 

c. 

...  SINGLE  INPUT 

0246 

c 

0247 

75 

UR I T  E ( IU, 1 40) 

0248 

READ  < IR,*)N1 , A I N  < 1 ) 

0249 

IF(N1.LT.1  .OR.  N1  .GT.  100)  URITEdU,  135) 

0250 

IF ( N1 . LT . 1  .OR.  N1  .GT.  100)  GO  TO  70 

0251 

D< ISYS,N1 )  =  A I N ( 1 ) 

0252 

GO  TO  70 

0253 

c 

0254 

c. 

...MULTI-INPUT-  10  NUMBER . 

0255 

c 

' 


0256 

0257 

0258 

0259 
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0261 

0262 

0263 

0264 

0265 

0266 

0267 

0268 

0269 

0270 

0271 

0272 

0273 

0274 

0275 

0276 

0277 

0278 

0279 

0280 

0281 

0282 

0283 

0284 

0285 

0286 

0287 

0288 

0289 

0290 

0291 

0292 

0293 

0294 

0295 

0296 

0297 

0298 

0299 

0300 

0301 

0302 

0303 

0304 

0305 

0306 
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30  UR  I T  E ( IU, 1 45) 

READdR,*)N1,(AIN(J),J=1,10) 

IF(N1 - G  T .90)  UR  I TE ( IU, 1 35) 

IF  ( N1  .GT.  90)  GO  TO  70 
DO  85  J-1 F 1 0 
JH1 =J+N1 -1 
D ( I  SYS , JH1 )=AIN(J) 

85  CONTINUE 
GO  TO  70 
C 

c 

90  CONTINUE 

C  CALL  CLOSEdDCB,  IERR) 

C  CALL  OPEN dDCB,  IERR, IF I LE, 0,77, 1 46, 1 70) 

CALL  RUNDF dDCB,  IERR) 

IF ( IERR  .GE.  0)  GO  TO  91 
STOP  4 

91  DO  95  J=1 ,6 
C 

c 

DO  820  I  =  1,40 
BUF(I)  =  D ( J,  I ) 

820  CONTINUE 

CALL  UR I TF (IDCB,IERR,BUF,80) 

DO  822  I  =  41,85 
IDUH  =1-40 
BUFdDUM)  =  D(J, I) 

822  CONTINUE 

CALL  URITFdBCB,  IERR  ,BUF, 90 ) 

I F ( IERR  .GE.  0)  GO  TO  97 
CALL  FHGER ( IERR) 

STOP  4 

97  CONTINUE 

U  R I T  E ( IU,305) J 
95  CONTINUE 

CALL  CLOSEdDCB, IERR) 

STOP 

100  FORMAT (5G12.5) 

105  FORMAT (17,1 0G1 2. 5) 

110  FORMAT ( 5G1 2. 5) 

115  FORMAT (5(G12.5, IX) ) 

120  FORMAT < 1 X , 60 ( 1 H-) ) 

125  FORMAT ( 1 X, ''MODEL  ,  SYSTEM  ,  CONTROL-ACTION  ,  RUN  ',/, 

1  IX/-  -  -  -) 

130  FORMAT ( 1 X/CHANGE  ,  CHANGES  ,  RESET  ,  <",2A1,I1/>  ",/, 

i  IX/-  -  -") 

135  FORMAT ( 1  X/URONG  NUMBER  OF  INDEX  TRY  AGIAN") 

140  FORMAT ( 1 X/ INPUT  INDEX  ,  VALUE") 

145  FORMAT ( 2X , "INPUT  INDEX  ,  VALUES  ....<><>") 

150  FORMAT ( "<01 >  TYPE  OF  SIMULATION  =" , 12 ) 

1  500  FORMAT  ( "  1 ",/////,"  ♦  *  *  " ,  A4/MQDEL  *  *  *',///) 


. 


184 


0307 

0308 

0309 

0310 

0311 

0312 

0313 

0314 

0315 

0316 

0317 

0318 

0319 

0320 

0321 

0322 

0323 

0324 

0325 

0326 

0327 

0328 

0329 

0330 

0331 

0332 

0333 

0334 

0335 

0336 

0337 

0338 

0339 

0340 

0341 

0342 

0343 

0344 

0345 

0346 

0347 

0348 

0349 

0350 

0351 

0352 

0353 

0354 

0355 

0356 

0357 


155  FORMAT ( 1 X, '<1 0>  N(  B1  )=',I2,/, 

1  ( 1 X , ' < 1 1 >  B1  =',5612. 5) 

160  FORMAT ( 1 X,'<20>  N(  B2  )=',I2,7, 


1  ( 1 X , '<21 > 

165  FORMAT (IX, '<30>  N( 
1  < 1 X,'<31> 

170  FORMAT ( 1X,'<40>  N( 
1  ( 1 X, '<41 > 

175  FORMAT ( 1 X, '<5Q>  N( 
1  ( 1 X, '<51 > 


B2 

D1 

D1 

Cl 

Cl 

A 

A 


=' ,5G1 2.5) 

)=',I2,/, 

=',5G12.5) 

)=',I2,/, 

=' ,5G12.5)  ) 

>  =  M2,/, 

=  ' ,5G12.5)  ) 

180  FORMAT ( 1 X,'<1 0>  T-RK  =',F2.0,'  FOR  G(11)',/, 

1  1X,'<11>  GAIN  =' ,G1 2.5, '<1 2-1 4>  T=',3G12.5) 

185  FORMAT ( 1 X, '<20>  T-RK  =',F2.0,'  FOR  G(12)',/, 

1  1 X , '<21 >  GAIN  =',G12.5,'<22-24>  T=',3G12.5) 

190  FORMAT ( 1  X, '<30>  T-RK  =',F2.0,'  FOR  G(  D>V, 

1  1 X , ' <3 1 >  GAIN  ='  ,G1 2.5, '<32-34>  T  =  ',3G12.5) 

195  FORMAT ( 1 X, '<40>  T-RK  =',F2.0,'  FOR  G(  C)',/, 

1  1 X , '<41 >  GAIN  ='  ,G1 2.5, '<42-44>  T=',3G12.5) 

200  FORMAT ( 1 X, '<02>  T-SA  =',G12.5, 

1  1 X, '<03>  RK-DT=',G12.5,'<04>  RK-NG=',F3.0,/ , 

2  1 X,'<05>  ST - D I = ' ,G1 2.5, '  <06>  N-LIM=' ,Gt 2.5) 

205  FORHAT ( 1 X, '<07>  K11=',F3.0,'  <08>  K12=',F3.0, 

3  1 X,'<09>  KD  =' ,F3.0, 1 X,F3.0) 

210  FORMAT ( '  <80>  ADC  Y=',G12.5,'  *  SADC  +  ',G12.5,/, 

S  '  <82>  DAC  U.S=',G12.5,'  ♦  UN  +  -,G12.5) 

2150  FORMAT ( ' 1 ',///// , '  *  *  *  ',A4, 'SYSTEM  *  *  *',///> 

215  FORMAT ( '<01>  NG1 =' , 12, '  <02>  NF  =',I2,'  <03>  ND  =',I2,/, 

1  1 X,  '<04>  NH  =' , 12 , '  <05>  NG2=',I2,/, 

•  2  1 X , ' < 0 6 >  K1  =' , 12, '  <07>  K2  =',I2,'  <08>  KD  s,I2,/, 

3  IX,'  FIX-N  =' , 12, '  <09>NDEL=',I2  ) 

220  FORMAT < 1 X, '<1 0>  G1  =',5G12.5) 

225  FORMAT ( 1 X, '<20>  F  =' ,5G1 2.5) 

230  FORMAT ( 1 X , '<30>  D  =',5G12.5) 

235  FORMAT < 1 X, '<40>  H  =' ,5612.5) 

240  FORMAT < IX,' <50>  G2  =' , 5G1 2.5) 

245  FORMAT ( 1 X, '<60>  TYPE  OF  IDENT=X, 12,'  ', A4,/, 

1  1 X,'<61>  R0U-1=',G12.5,'  <62>  R0U-2=' ,G1 2.5,/, 

2  1X,'<63>  P-INT  =  ' ,G12.5, '  <64>  I -F I X  = ' , 2 1 3 , / , 

3  1X,'<66>  G*COP=' ,G1 2.5, '  <67>  L-PHI=' ,2G1 2.5) 

250  FORMAT ( 1 X, '<70>  TYPE  OF  SET-P0INT=',F2.0,/, 

1  1 X, '<71 >  HIEGHT  =',G12.5,'  <72>  LEINGHT=',G1 2.5,/, 

2  1X,'<73>  START  =' ,G1 2.5, '  <74>  FINAL  =',G12.5,/, 

3  1X,'<75>  BASE  =' ,G1 2.5) 

255  FORMAT (IX,  '<80>  TYPE  OF  DISTURBANCES, F2.0,/, 

1  1 X,'<81>  HIEGHT  =',G12.5,'  <82>  LEINGHT=X,G12.5,/, 

2  1X,'<83>  START  =',G12.5,'  <84>  FINAL  =',G12.5,/, 

3  1X,'<85>  BASE  =',G12.5) 

2600  FORMAT ( ' 1 ',/////,'  *  *  ♦  ' ,A4, 'CONTROL  *  *  *',///) 

260  FORMAT ( '<01 >  TYPE  OF  CONTROLLERS , I2,4X,A4,/, 

1  1 X, '<02>  LOUER-U=' ,G1 2.5, '  <03>  UPPER-U=' ,G1 2.5,/, 

2  1 X, '<04>  PERCENT  =  ' , G 1 2.5, '  USTEP  =',G12,5) 


< 


. 


0358  265  FORMAT (IX, x < 0 5 >  TYPE  OF  PID=" , 12, 1  OX/  0  — >  Z-TRAN",/, 

0359  1  IX,  28X,  ■'  1  — >  S-TRAN",/, 

0360  2  1  X,  "<06>  PID-ACTIQN="  ,G1 2.5) 

0361  270  FORMAT ( 1  X, "<07>  KP  =",G12.5,/, 

0362  1  1 X,  "<08>  KI  =  "  ,G1 2.5,/, 

0363  2  1  X ,  "<09>  KD  =",012.5) 

0364  275  FORMAT ( 1 X,  "Q-FILTER _ " ,/, IX, '<1 0>  NQ  =",I2,/, 

0365  1  1X,"<11>  QD(Z)=",5G12.5,(/,10X,5G12.5>) 

0366  280  FORMAT (1X,"<15>  QN(Z)=" ,5G1 2.5, (/, 1 0X,5G1 2.5) ) 

0367  285  FORMAT ( 1 X, "P-FILTER _ ",/, 

0368  X  1 X , "<2Q>  TYPE  =",I2,10X,"  0  — >  Z-TRAN",/, 

0369  X  IX,"  T-S  ="  ,G1 2.5  ,"  1  -  — >  S-TRAN") 

0370  290  FORMAT ( 1 X , " <2 1 >  T-RK=",I2,/, 

0371  X  1 X ,  "<22>  GAIN  =",G12.5,"  <23-25>  T=",3G12.5) 

0372  295  FORMAT ( 1 X, "<3Q>  NP  =",I2,/, 

0373  i  1 X, "<31 >  PN  =", 4G12.5,/, 

0374  X  1  X, "<35>  PD  =" , 4G1 2.5) 

0375  300  FORMAT ( 1 X, "R  FILTER _ ",/,1X,"<50>  T-RK=" ,F2.0,/, 

0376  1  1 X, "<51 >  GAIN=" , G1 2.5, "  <52-53>  T=" ,2G1 2.5,"<54>  K=",F2. 

0377  301  FORMAT ( 1 X, "<60>  U-STEADY-STATE=" ,G1 2.5) 

0378  302  FORMAT (IX,"  PID  SETTINGS  ( AKPS, AKIS, AKDS) ",/, 

0379  X  IX,"  POSITIONAL  OR  INCREMENTAL  CONTROLLER",/, 

0380  1  1 X, "<61 -63>  " , 3 G 1 2.5) 

0381  303  FORMAT (IX,"  Q-FILTER  SETTINGS  (AKP,AKI,AKD)" ,/, 

0382  X  1  X, "<64-66>  "  ,3G1 2.5) 

0383  305  FORMAT ( "  SECTOR  " ,  12 , "  HAS  BEEN  URITTEN  ") 

0384  310  FORMAT < 2A1 ,A2) 

0385  315  FORMAT (II) 

0386  END 


* 

;« 
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SERSUM  T=00004  IS  ON  CR00146  USING  00011  BLKS  R=0087 
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O04B 


FTN4 

PROGRAM  ERSUM (  ,101) 

C 

DIMENSION  IDCB(144),BUF(56),ISTR(20) , 

*  II (20), 12(20), IR(20), 

*  IFILE(3) ,LU(5) , 

*  FEED1 (20) ,FEED2<20) ,SIAE(20) 

C 

CALL  RMPAR(LU) 

LUN=LU  <  1 ) 

C 

C 

URITE(LUN, 1 000) 

READ(LUN, 1 002)  (IFILE(I) , 1=1 ,3), ICR 
C 

URITE(LUN, 1 006) 

READ(LUN,*)  MAXREG 
C 

URITE ( LUN , 1010) 

READ(LUN,*)  ILL) 

C 

C 

BO  50  IREG  =  1  , MAXREG 
URI TE < LUN ,1012)  IREG 
READ (LUN,*)  II (IREG), I2( IREG), IR (IREG) 
50  CONTINUE 
C 

c 

c 

c 

CALL  OPEN(IDCB,IERR,IFILE, 0,77, ICR, 144) 
IF(IERR.GE.O)  GOTO  70 
CALL  FMGER ( I ERR ) 

STOP 

C 

C-  READ  HEADER  RECORD. 

C 

70  CONTINUE 

CALL  READF( IDCB, IERR,BUF,40, ILEN) 

IF (IERR.LT.O)  CALL  FMGER ( I  ERR ) 

UR I T E ( ILU, 1 020)  (BUF(I) ,1=1 ,20) 

DO  900  IREG  =  1 , MAXREG 
C 

IBEG  =  II (IREG) 

IEND  =  0 
I DUM  =  IREG  -  1 

IF(IREG  .NE.  1)  IEND  =  ISTR(IDUH)  +  IR(IDUM) 
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■ 
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IRU  =  IBEG  -  IEND  -1 
CALL  POSNT ( IDCB, IERR, IRU,0) 

CALL  HSI AE ( IDCB, BUF , ISTR* IREG ) , BEG, END , 

*  SI AE ( IREG), FEED 1( IREG) ,FEED2( IREG), 

*  II <  IREG) ,I2(IREG) ,IR( IREG)) 

URITE ( ILU, 1  022 )  IREG,  BEG , END 

1022  FORMAT * *  REGION  ",14,'  BEGINS  AT  ",F 7.2/  ENDS  AT  ',F7.2> 

C 

C 

900  CONTINUE 
C 

C-  PRINT  OUT  SI AE  RESULTS. 

C 

CALL  RUNDF ( IDCB, IERR) 

CALL  READF( IDCB, IERR, BUF, 40, ILEN) 

U R I T E ( I LU , 95 1 )  (BUF(I) ,1  =  1 ,20) 

C 

DO  950  IREG  =  1 , HAXREG 
U R I T E ( I LU , 953 )  IREG, 1 1 ( IREG ) , 12 ( IREG) , 

*  ISTR( IREG) ,FEED1 (IREG), FEED2( IREG), 

*  IR ( IREG ) , SI AE ( IREG ) 

950  CONTINUE 
C 

CALL  CLOSE( IDCB, IERR) 

STOP 

C 

c - 

c- 

C-  FORMAT  STATEMENTS. 

C- 

c 

951  FORMAT (///  *  *  FOR  FILE  ON  *  *  ",/,20A4,/) 

953  FORMAT*//  FOR  REGION  ",I2,5X/  FROM:  ',14/  TO:  ',14, 

*  //  FEED  CHANGE  ATSH/X/  FROM:  SF6.3,"  TO:  \F6.3 

*  //  SIAE  AFTER  ",12/  ITERATION  IS  =",F10.4) 

1  000  FORMAT*//  ENTER  FILENAME ( 3A2 )  AND  CARTRIDGE* 13) :  '  / ) 

1002  FORMAT ( 3A2 , 13 ) 

1006  FORMAT*//  ENTER  NO.  OF  REGION  <MAX  =  20):  "/_•) 

1010  FORMAT*//  ENTER  OUTPUT  DEVICE:  "//) 

1012  FORMAT*//  FOR  REGION:  ",I2, 

*  //  ENTER  START, FINAL  AND  NO.  OF  ITERATIONS:  V'/) 
1020  FORMAT*//  PLOTTING  FOR  ",/,20A4) 

END 


SHPL01  T=00004  IS  ON  CR00136  USING  00028  BLKS  R=0000 
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C 

c*  THIS  PROGRAM  PLOTS  DATA  STORED  IN  DATA  FILE. 

C 

FTN4 

PROGRAM  HPLOT (  ,101) 

C 

DIMENSION  IDCB<144),BUF(56),X(?02),Y(702), 

*  11(15), 12(15), IR(15), 

*  SPAN(6),ZER0(6), IFL0C6) ,Y00(6) , 

*  IP0S(11),SSTART(11,15), RRANGE ( 1 1,15) 

*  NPAR(5),IFILE(3),LU(5), 

*  DSTART (6) ,DRANGE(6) 

C 

DATA  IPOS/9, 10, 12, 11 ,13,14,17,0,0,0,0/, 

*  CON A/ 1 .0/,C0NB/0./, 

*  IFLO/O, 1,1, 0,1,1/, 

*  SPAN/. 12857, 2. 2975, 2. 551 7,1 .,2.50,1 .3140/ 

*  ZER0/89.167,0. ,0.  ,0.  ,0.  ,0./, 

*  SSTART/1 10*0./, RRANGE/1 10*0./, 

*  Y00/0.,4.,5.5,0.,4.,5.5/, 

*  DSTART/93. 5, 0.00, 12. 0,1 .00,0.00,0.00/, 

*  ORANGE/ 1 .00,0.00,12.00,3.00,0.00,0.00/ 
CALL  RMPAR(LU) 

LUN=LU ( 1  ) 

C 

C*  OBTAIN  INFORMATION  ON  WHERE  DATA  IS  STORED 
C 

WRITE ( LUN ,  1 000 ) 

READ(LUN, 1 002)  (IFILE(I) ,1=1 ,3) ,ICR 
C 

WRITE ( LUN , 1 004 ) 

READ(LUN,*)  ( NPAR ( I ) , 1=1 ,5) 

C 

WRITE ( LUN , 1 006 ) 

READ ( LUN , * )  MAXREG 
C 

WRITE (LUN , 1 008 ) 

READ(LUN,*)  IYOPT 
C 

WRITE(LUN, 1010) 

READtLUN,*)  ILU 
C 

WRITE (LUN, 101 1 ) 

READ ( LUN , ♦ )  ST 

C 

C*  READ  IN  INFORMATION  FOR  EACH  REGION. 

C 


•  g 


. 

*  '.a  .  j  : 
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DO  50  IREG  =  1  , MAXREG 

URITEdUN, 1012)  IREG 

READ (LUN,*)  II (IREG) , 12 ( IREG) fIR( IREG) 

C 

URITEdUN, 1014) 

READdUN,*)  IDUM 
IF ( IDUM.NE. 1 )  GO  TO  20 
C 

C*  OBTAIN  RANGE  AND  START  PONT  OF  EACH  REGION  FOR  PROCESS  MEASURE 
C 

URITE(LUN,  1016)  MAXREG 
READdUN,*)  (SSTARTd, IREG), 1  =  1 ,6) 

URITE(LUN, 1 017)  MAXREG 
READdUN,*)  (RRANGE( I , IREG) , 1  =  1 ,6) 

C 

20  IF ( IDUM.NE. 2)  GO  TO  22 

C 

C-  TAKE  DEFAULT  SCALING. 

C 

DO  24  1=  1,6 

SSTART ( I , IREG)  =  DSTART (I ) 

RRANGE ( I ,  IREG )  =  DRANGE(I) 

24  CONTINUE 

C 

22  URITE(LUN, 1018) 

READdUN,*)  IBUH 
IF < IDUM.NE. 1 )  GO  TO  30 
D 

C*  UANT  START  AND  RANGE  OF  PARAMETERS. 

C 

URITE ( LUN, 1016)  MAXREG 
READdUN,*)  (SSTART(I,IREG),I=7,1 1 ) 

URITE (LUN, 1017)  MAXREG 

READ (LUN,*)  (RRANGE (I, IREG), 1=7, 11  ) 

C 

30  CONTINUE 

C 

50  CONTINUE 
C 

C*  OBTAIN  THE  STARTING  POSITION  OF  EACH  TYPE  OF  PARAMETER 
C 

DO  60  I  =  8,11 
IDUM  =1-1 
IDU1  =1-7 

IPQS(I)  =  IPOS ( IDUM )  +  NPAR(IDUI) 

60  CONTINUE 
C 

c 


CALL  OPEN ( IDCB, IERR, IFILE, 0,77, ICR,  144) 
IF (IERR.GE.O)  GOTO  70 
CALL  FMGERdUN,  IERR) 


* 


. 

i 
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0100 

0101 
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0116 
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0120 

0121 

0122 

0123 

0124 


STOP 

C 

C-  READ  HEADER  RECORD. 

C 

70  CONTINUE 

CALL  READF(IDCB,IERR,BUF,40,ILEN) 

IFUERR.LT. 0)  CALL  FHGER(LUN, IERR) 

URITE(LUN, 1 020)  (BUF(I) ,1=1 ,20) 

DO  900  IREG  =  1 , HAXREG 
C 

URITE(LUN,2040) 

IBEG  =  II ( IREG) 

IEND  =  0 

IF  (IREG  .NE.  DIEND  =  I2UREG-1) 

IRU  =  IBEG  -  IEND  -  2 
CALL  P0SNT(IDCBfIERR,IRU,0> 

CALL  LOCF(IBCB,IERR,IREC,IRB,IOFF> 

C 

C*  FILL  IN  TIME  INFORMATION  (X-AXIS). 

C 

J  =  0 

DO  100  I  =  II (IREG) ,I2(IRE6>,  IR(IREG) 
J=  J  +  1 

X(J)  =  FLOAT ( I - 1 1 (IREG))*ST 
100  CONTINUE 
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C 

C*  COHPUTE  NO.  OF  POINTS  TO  BE  PLOTTED 
C 

NPT  =  I2(IREG)  -  1 1 (IREG)  *  1 
C 

C*  INITIALIZE  PLOTTING. 

C 

CALL  PLOTS ( 0 - 0 , ILU) 

IF ( ILU.EQ.4)  CALL  FACTOR*. 6) 

CALL  PLOT (1.87,1.61,-3) 

CALL  HCALE ( X , 5 . 0 , NPT , 1 ) 

C 

C*  DETERMINE  WHICH  OPTION  TO  PLOT. 

C 

IBEG  =  1 
IEND  =  3 

IF ( IYOPT.EQ. 1  -OR.  IY0PT.EQ.3)  GOTO  150 
C 

IBEG  =  4 
IEND  =  6 

IF (IYOPT.EQ. 2)  GOTO  150 
C 

IBEG  =  1 
150  CONTINUE 

C 

c 
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C-  START  PLOTTING  THE  FLOU  MEASUREMENT. 

C 

DO  200  I  =  I  BEG  r I  END 
ITYP  =  I 
C 

c-  IF  ONLY  PLOTTING  BOTTOM,  PLOT  FEED  INSTEAD  OF  REFLUX. 

C 

IF ( IY0PT.EQ.2  .AND.  ITYP.EQ.6)  ITYP  =  3 
IF (IY0PT.EQ.3  .AND.  ITYP.EQ.2)  ITYP  =  4 
START  =  SSTART (ITYP,IREG) 

RANGE  =  RRANGE( ITYP, IREG) 

C 

C*  CALL  SUBROUTINE  TO  PLOT  OUT  MEASUREMENT  FOR  EACH  REGION. 

C 

CALL  HPLT(IDCB,BUF,X,Y,IREC,IRB,IOFF, 

*  II (IREG), 12 (IREG), IR< IREG), BEG,  END, 

*  IPOS < ITYP ), SPAN (ITYP), ZERO ( ITYP >,YOO< I) 

*  IFLO( ITYP), START, RANGE, ITYP) 
IFUTYP.EQ.3.AND.IEND.EQ.6.AND.IY0PT.EQ.4)  CALL  PLOT(0.,-9 

URI TE ( LUN , 2000 )  ITYP, IREG 
URITE(LUN,2020)  BEG, END 

2020  FORMAT ( '  PLOTTING  FROM  ',F7.2/  TO  ",F7.2) 

200  CONTINUE 

C 

C*  PLOT  OUT  PARAMETERS  FOR  EACH  REGION. 

C 

IF < I LU  .EQ.  4)  GO  TO  800 
CALL  PLOT (0. ,-9.5, -3) 

CALL  PLOT ( 1 .397, .2733,-3) 

C 

C 

DO  300  I  =  1,5 
ITYP  =6+1 


C 

START  =  SSTART ( ITYP, IREG) 

RANGE  =  RRANGE< ITYP, IREG) 

C 

C*  CALL  SUBROUTINE  TO  PLOT  OUT  PARAMETERS  FOR  EACH  REGION. 

C 

CALL  HPLT(IDCB,BUF,X,Y,IREC,IRB,IOFF, 

*  II (IREG ),I2< IREG), IR( IREG), BEG, END, 

*  IPQS(ITYP),CQNA,C0NB,Y00(1),NPAR(I),O, 

*  START, RANGE, ITYP) 

URITE(LUN, 2000)  ITYP,  IREG 
URITE(LUN,2020)BEG,END 

C 

300  CONTINUE 

CALL  PLOT (-1.397, -.2733, -3) 

800  CALL  PLOT (0. ,0. ,999 ) 

C 

c 


. 
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900  CONTINUE 


CALL  CLOSE*IDCB,IERR) 
STOP 


C- 

C-  FORHAT  STATEMENTS. 

C- 

C 

1000  FORMAT (//  ENTER 
1002  FORMAT < 3A2, 13 ) 
FORMAT*//  ENTER 
FORMAT*// 
FORMAT*// 
FORMAT*// 
FORMAT*// 


FILENAME(3A2)  AND  CARTRIDGE ( 13) 


t  ~ 


') 


1004 

1006 

1008 

1010 

1011 


) 


ENTER 

ENTER 

ENTER 

ENTER 


NO.  OF  G1,F,D,H,G2  PARAMENTERS:  J 
NO.  OF  REGION  < MAX=  1 0 ) :  "//) 

1  =  TOP,  2  =  BOTTOM,  3  =  FD,TC,BC,  4  =  ALL: 
PLOT  DEVICE  ( 4=HP2648  6=LP):  "/  "> 

SAMPLE  TIME  IN  MINUTES:  ") 

1012  FORMAT*//  FOR  REGION:  ',12, 

*,//  ENTER  START,  FINAL  AND  REPEAT  FACTOR:  "/  "> 

1014  FORMAT*//  ENTER:  O-AUTOSCALE  1-MANUAL  ENTRY  2-DEFAULT:  "/ 

1016  FORMAT ( "  ENTER  ",12/  START  VALUE  FOR  REGIONS,  0  =  AUTO  SCAL 

1017  FORMAT*"  ENTER  " ,  12,"  RANGE  FOR  DIFFERENT  REGION  (UNITS/TICK 

1018  FORMAT*//  UANT  TO  ENTER  SCALE  FACTOR  FOR  PARAM.", 

*  "(1  -  YES):  _") 

1020  FORMAT*"  *  *  PLOTTING  FOR  *  *",/,20A4) 

2000  FORMAT*//  *  *  COMPLETED  TYPE  ",12/  FOR  REGION  ',12/  *  *" 
2040  FORMAT*/) 

END 
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0001 

C 

0002 

C* 

THIS  SUBROUTINE  DOES  THE  ACTUAL  PLOTTING  FOR  EACH 

0003 

C 

0004 

FTN4 

0005 

SUBROUTINE  HPLT < IBCB,BUF , X , Y , IREC, IRB f IOFF , 

0006 

*  11 ,12, IR  ,BEG,END, 

0007 

*  IPOS, SPAN, ZERO, YOO,NPAR,IFLG, 

0008 

*  START, RANGE, ITYP) 

0009 

C 

0010 

DIMENSION  IBCB( 1 ) ,BUF ( 1 ) ,X( 1 ) ,Y( 1 ) 

0011 

DATA  ILEN/1 1 2/, IFLG/1 / 

0012 

IF ( NPAR  .EQ.  0)  GO  TO  500 

0013 

C 

0014 

C 

0015 

NPT  =12-11+1 

0016 

K  =  0 

0017 

C 

0018 

C 

0019 

DO  100  II  =  1 ,NPAR 

0020 

J  =  0 

0021 

CALL  APOSN(IDCB,IER,IRCD, IRB, IOFF) 

0022 

CALL  READF ( IDCB,IERR,BUF,20,ILEN) 

0023 

I  UAL  =  IPOS  +  K 

0024 

K  =  K  +1 

0025 

C 

0026 

C 

0027 

DO  110  I  =  1 ,NPT 

0028 

CALL  READF (IDCB , IERR,BUF, 112, ILEN ,LEN) 

0029 

IF ( ILEN  .EQ.  -1)  GO  TO  900 

0030 

IF ( I  .EQ.  1)  BEG  =  BUF (3) 

0031 

IF < I  .EQ.  NPT)  END  =  BUF(3) 

0032 

C 

0033 

J  =  J  +  1 

0034 

Y(J)  =  BUF ( I UAL ) 

0035 

IF ( IFLO.EQ. 1 )  Y ( J)  =  SQRT(Y(J>> 

0036 

Y ( J )  =  SPAN  ♦  Y ( J )  +  ZERO 

0037 

110 

CONTINUE 

0038 

C 

0039 

C 

0040 

IF (K  .NE.  4)  GO  TO  120 

0041 

CALL  PLOT ( 0. , 0 . , 1 0) 

0042 

CALL  PLOT (0. — 3) 

0043 

c 

0044 

120 

CONTINUE 

0045 

X0=0 

0046 

IF ( YOO  .EQ.  0.)  GO  TO  200 

0047 

YO  =  YOO 

0048 

GOTO  210 

'  -  3 
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C 


200 

YO  =  3.5 
IFdI.EQ.I  .OR. 

II.EQ.4)  YO  =  0 

210 

CALL  PLOT ( XO, YO, 

-3) 

C-  START  PLOTTING. 


C 


C 


C 


C 


C 


C 


C 


C 


C 


N=NPT 

CALL  SYMBOLd .88, -.55, .15, 10HTIME  (MIN), 0.0, 10) 
CALL  AXIS(0.,0.,20H  , 

*  -20,5.0,0.0,X(N+1),X(N+2)) 

Y (N  +  1 )  =  START 
Y (N  +  2)  =  RANGE 

IF ( START  .EQ.  0.)  CALL  HCALE ( Y,3.0,NPT, 1 ) 
GOTOd ,2,3,4,5,6,7,8,9,10,11 )  ITYP 


1  CALL  SYMBOLd. 4, .975,. 15, 10HXB  -  LUTZ], 90.0, 10) 
CALL  AXIS(0.,0.,0,0,3.0,90.0,Y(N+1),Y(N+2)) 

GOTO  250 

2  CALL  SYMBOLd. 4, .975,.  15, 1  0HRE  -  CG/S] ,90.0,10) 
CALL  AXIS(0.,0.,0,0,3.0,90.,Y(N+1),Y(N+2)) 

GOTO  250 


3  CALL  SYMBOL (-.4, .975, .15,1 OHFE  -  CG/S3 , 90.0,10) 
CALL  AXIS(0.,0.,0,0,2.0,90.,Y(N+1),Y(N+2) ) 

GOTO  250 

4  CALL  SYhB0L<-.4, .975, . 1 5, 1 OHXB  -  CUT%3 ,90-0,10) 
CALL  AXIS(0.,0.,0,0,3.0,90.,Y(N+1),Y(N+2)) 

GOTO  250 

5  CALL  SYMBOLd. 4, .975,.  1  5, 1  0HST  -  CG/S] ,90.0, 1  0 ) 
CALL  AXIS(0.,0.,0,0,3.0,90.,Y(N+1),Y(N+2)) 

GOTO  250 


6  CALL  SYMB0L(-.4,.975,.15,10HDI  -  CG/SD , 90.0,10) 
CALL  AXIS(0.,0.,0,0,3.0,90.0,Y(N+1),Y(N+2) ) 

GOTO  250 

7  CALL  SYMBOLd. 4, .975,.  15, 10HPARAM  -  G1 , 90. 0,10) 

CALL  AXIS ( 0 .  ,0 . , 1 8H  ,18,3.0,90., 

*  Y(N+1),Y(N+2)> 

GOTO  250 


8  CALL  SYMBOLd. 4, .975,.  15, 10HPARAH  -  F  ,90.0,10) 
CALL  A X I S < 0 . ,0. , 1 8H  ,18,3.0,90., 

♦  Y (N+1 ) ,Y (N+2) ) 

GOTO  250 


. 

l 

r 
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0115 
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0123 
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9 

CALL 

SYMBOL (-.4, .975,. 15,1 OHPARAM 

-  D  ,90.0,10) 

CALL 

AXIS ( 0.  ,0. , 1 8H 

,18,3.0,90., 

♦ 

Y (N+1 ) , Y (N+2) ) 

GOTO 

250 

10 

CALL 

SYMBOL*-. 4,. 975,. 15,1 OHPARAM 

-  H  ,90.0,10) 

CALL 

A X I S ( 0 .  ,0. , 18H 

,18,3.0,90., 

♦ 

Y (N+1 ) , Y (N+2) ) 

GOTO 

250 

11 

CALL 

SYHBQL(-.4,.975,.15,1 OHPARAM 

-  G2, 90. 0,10) 

CALL 

A  X I S ( 0 . ,0. , 1 8H 

,18,3.0,90., 

* 

Y(N*1 ) ,Y (N+2) ) 

C 

250  CONTINUE 

CALL  LINE ( X , Y , NPT ,1,0,1) 

100  CONTINUE 
C 
C 

IF ( ITYP.EQ.3  .OR.  ITYP.GE.6)  CALL  PLQT( 0 . , 0 . , 1 0 ) 
C 

C-  CHECK  IF  UE  HAVE  TO  RESET  THE  ORIGIN. 

C 

IF (NPAR.EQ.  1  .OR.  NPAR.EQ.4)  YO  =  0 
IF (NPAR.EQ.2  .OR.  NPAR.EQ. 5)  YO  =  -3.50 
IF (NPAR.EQ. 3  .OR.  NPAR.EQ. 6)  YO  =  -7. 

CALL  PLOT (0. , Y0,-3) 

C 

500  CONTINUE 
RETURN 
C 
C 

900  URITE(LUN,999) 

999  FORMAT ( '  REACHED  END  OF  FILE  ') 

CALL  CLOSE(  IDCB, IERR) 

STOP 

END 
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FTN4 

PROGRAM  GCLNK  (3,70),  GC  DATA  TRANSFER,  790301,  VFB 
C 
C 

C  THIS  PROGRAM  IS  USED  TO  RECEIVE  GC  DATA  FROM  A  REMOTE 

C  CPU.  IT  USES  CLASS  I/O  TO  RECEIVE  THE  DATA. 

C  THE  PROGRAM  IS  ALUAYS  READY  TO  RECEIVE  THE  DATA. 

C  THE  PROGRAM  STORES  DATA  IN  REAL  TIME  COMMON  AREA. 

C  START  OF  THE  AREA  MUST  BE  SET  UP  IN  REAL  TIME 

C  COMMON  POINTER  AREA.  CURRENTLY  POINTER  4,  WORD  8 

C  OF  RT  COMMON  IS  USED  TO  POINT  TO  THE  GC  DATA  AREA. 

C  THE  LENGTH  OF  THE  AREA  IS  80  UORDS. 

C 

C 

c 

c 

COHMON  /  /  ICOM(I) 

DIMENSION  IBUF (50,10),  ICLAS(IO) 

C 

C  INITIALIZE  I/O,  CHECK  STATUS,  LOCK  RECEIVING  UNIT, 

C  START  CLASS  I/O  AND 

C  SET  UP  COMMON  ADDRESSING 

C 

DATA  LUGC/33/ ,  NBUF/10/,  NEXT/O/ 

DATA  IRLEN/50/ ,LEN/50/ 

DATA  ICLAS/1 0*1 OOOOOB/ 

IAS=ICQM ( 8 ) 

IAD=IAS+6 
ICQM( IAS)=0 
LU=LOGLU(ISESN) 

CALL  EXEC  ( 1 3,LUGC, IST1,IST2,IST3) 

IF  (IST3.LT. 0)  GO  TO  800 
CALL  L UR Q ( 1  ,LUGC,1 ) 

C  CALL  EXEC  (22,1) 

DO  100  1=1 ,NBUF 
1 1  =  I 

CALL  EXEC ( 17,LUGC, IBUF( 1,1)  , IRLEN , I  PI , IP2,ICLAS( I) ) 
CALL  ABREG  ( I  A, IB ) 

D  WRITE  ( LU , 1 030 )  IA, IB 

D1030  FORMAT  ("  ALOCATE  1:  IA,IB=",208) 

IF  ( IA.LT.O)  GO  TO  110 
I CL AS  < I )  =  IAND( ICLAS( I) ,07777?B) 

100  CONTINUE 
C 

110  NBUF=II-1 

WRITE (LU, 1 035)  NBUF 

1035  FORMATdX, "NUMBER  OF  BUFFERS2" ,13) 

C  ICLAS=IOR  ( ICLAS ,20000B ) 


0049 

0050 

0051 

0052 

0053 

0054 

0055 

0056 

0057 

0058 

0059 

0060 

0061 

0062 

0063 

0064 

0065 

0066 

0067 

0068 

0069 

0070 

0071 

0072 

0073 

0074 

0075 

0076 

0077 

0078 

0079 

0080 

0081 

0082 

0083 

0084 

0085 

0086 

0087 

0088 

0089 

0090 

0091 

0092 

0093 

0094 


197 


C  I  CLAS  =  I  AND ( ICLAS ,077777B) 

C 

c 

C  START  OF  INPUT  LOOP 

C 

c 

500  CONTINUE 
C 

DO  200  IBUFN=1,NBUF 

CALL  EXEC(21 ,  ICLAS (IBUFN) ,IBUF( 1 , IBUFN ),LEN, IP  1 ,IP2,IP3) 

CALL  ABREG< IA, IBUFL) 

D  WRITE  <LU,  1 050)  IA, IBUFL 
D1050  FORMAT  ( "  GET  Is  IA,  IBUFL  =  '',208) 

I CLAS  < IBUFN)=0 
IF  ( IBUFL.EQ.O)  GO  TO  130 
IF  ( IBUFL. LT.O. OR. IBUFL. GT. 40)  GO  TO  150 
D  URITE  (LU,1060)  ( IBUF ( II I , IBUFN) , 1 11  =  1 , IBUFL) 

D1060  FORMAT ( 1 X,20A2) 

IF  (IAD+IBUFL.GT.IAS+IC0M(9)-1)  IAD=IAS+6 
D  URITE (LU, 1070)  LUGC, IRLEN,LEN, IBUFN, IAD 
1070  FORMAT ( "LUGC, IRLEN,LEN, IBUFN, IAD="/ 

25110) 

C 

DO  120  J=1 , IBUFL 
I  COM ( IAD  +  J)  =  IBUF( J, IBUFN) 

120  CONTINUE 
C 

130  IF  (ICOM(IAS).EQ.O)  GO  TO  140 
NEXT =NEXT  + 1 

IF  (NEXT.GE.NBUF)  GO  TO  999 
GO  TO  200 
C 

140  I AD=I AD+IBUFL 

150  CALL  EXEC  ( 1 7, LUGC , IBUFd , IBUFN ), IRLEN, IP1 , IP2 , ICLAS ( IBUFN ) ) 
CALL  ABREG  <  IA, IB) 

D  URITE  (LU, 1 030)  IA, IB 
200  CONTINUE 
GO  TO  500 
C 

C  END  OF  INPUT  LOOP 
C 

c 

800  CALL  ABORT  (0, 1 ,25H***GC  LINK  UNIT  DOUN  ***> 

999  CALL  LURQ(0,LUGC, 1 ) 

STOP 

END 
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FTN4 

PROGRAM  GCSCH(3,80) ,  PROGRAM  TO  INITIATE  GC  OH  COLUMN  HYL790 
COMMON  I  COM ( 1 ) 

DIMENSION  INAME1 (3) , I  NAME  2  <  3 ) 

DATA  INAME1 /2HGC,2HSU, 2HT  /, INAME2/2HSE,2HGS,2HT  / 
ISTR*IC0M(8) 

ILEN=IC0M(9)-7 
LUN=ICOM( ISTR+3) 

IHAX=ICQM( ISTR+4) 

IDFLT  =  I  COM ( I STR  +  5) 

C 

CALL  BOTCMdSTR,  ILEN,LUN, IHAX ) 

CALL  DQL(1,1,4,4,IERR> 

IF( IERR.NE. 1 )  GO  TO  500 
C 

C*  SCHEDULE  GCSUT  TO  RESET  GC  START  BUTTON. 

C 

CALL  EXECd  0  +  1 OOOOOB,  INAHE1 ) 

GO  TO  950 
STOP 
C 

C*  IF  ERROR  DETECTED,  TURN  PROGAM  OFF  AND  SET  COMPOSITION  TO  DEFA 
C 

500  URITE(LUN,600) IERR 

600  FORMAT ( "PROGRAM  GCSCH  FAILED  TO  CLOSE  DITIGAL  SUITCH,  IERR  - 
GO  TO  999 
C 

950  URITE(LUN,952) 

952  FORMATCSCHEDULE  OF  GCSUT  UNSUCCESSFUL" ) 

C 

999  IC0H(ISTR+1)  =  IDFLT 
CALL  EXEC< 1 0, INAHE2) 

STOP 

END 

C 

C 

c 

SUBROUTINE  BOTCMdSTR,  ILEN,LUN,  INAX ) 

COMMON  I  COM ( 1 ) 

DIMENSION  C0N(6) ,  IVAL(4),  ITIME (5) 

DATA  CON/IOO. ,1.,.1,. 001, 10.,. 01/, IBLANK/2H  / 

ICHR=2H  6 
C 

C*  SEARCH  FOR  CHARACTER  STRING  "4".  THE  COMPOSITION  OF  UATER  UIL 
C*  PRECEDE  THIS  STRING. 

C 

IADD-ISTR+6 
DO  10  1=1 , ILEN 
IPRE=I 

IFCICOMdADD  +  IPRE) .EQ. ICHR)  GO  TO  12 


. 
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10  CONTINUE 
60  TO  500 
C 

C*  STORE  ASCII  REPRESENTATION  OF  COMPOSITION  INTO  ARRAY  IVAL. 

C 

12  DO  20  1=1,4 
IPOS=IPRE-I 

IF ( IPOS.LT .0 )  IP0S=ILEN+IPQS+1 
IVAL(5-I)=IC0H(IADD+IPQS) 

20  CONTINUE 
C 

C*  CONVERT  THE  VALUE  IN  ASCII  INTO  A  DECIMAL  NUMBER. 

C 

VALU=0. 

DO  30  1=1,4 

IDUM=IAND( ISHFT ( IVAL( I ) ,-8) ,377B) 

IF( IDUM.EQ.32)  C0N(I)=0. 

VALU=VALU+CON( I )*FLOAT ( I DUN-48) 

IF( (I.EQ.2)  .OR. < I.EQ.4) )  60  TO  30 

IDUH=( 1-1 >/2+5 

IDU1 =IAND< IVAL( I ) ,377B) 

I F ( I  DU 1 .EQ.32)  CON( IDUM)=0. 

VALU=VALU+CON( I DUH ) ♦FLOAT ( I DU 1-4 8) 

30  CONTINUE 
C 

c 

C*  CONVERT  THE  COMPOSITION  TO  AN  INTEGER  BETUEEN  0-32767  FOR  METH 
C*  NOTE:  1Z=1000. 

C 

VALU=1 OO.-VALU 
IHET=IFIX(VALU*1000.) 

IF  (LUN.NE.34)  URITEtLUN, 1 03)  VALU 
103  FQRHAT  ( 1 OX,'*XB  =  M,F6.3) 

C 

C*  CHECK  IF  THE  VALUE  HAS  CHANGED. 

C 

IF< IMET.EQ.ICOH( ISTR+1 ))  GO  TO  200 
IF ( I  MET  .LT.  0)  GO  TO  220 
C 

c-  CHECK  THAT  GC  READING  HAS  NOT  CHANGE  BY  MAX  MAX  IS  STORED  IN 
C-  UORD  OF  GC  AREA. 

C 

IDIFF  =  I ABS( I COM ( I STR  +  1)  -  IMET) 

IFUDIFF.  GT.  ICOMUSTR  +  6)  )  GOTO  300 
C 

I  COM ( I STR  +  2)  =  0 
ICQM( ISTR+1 )=IMET 
C 

C-  BLANK  OUT  GC  AREA  IN  RT  COMHON. 

C 


DO  25  1=1 , ILEN 
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200 


II  =*  I 

ICOMUADD  ♦  II)  =  IBLANK 
25  CONTINUE 
C 

c 

RETURN 

C 

C*  IF  GC  READING  REMAINS  THE  SAME  FOR  IMAX  TIME,  GC  COMPUTER  IS  P 
C*  DOWN,  PUT  OUTPUT  SEGMENT  OF  STEAM  FLOW  ON  MANUAL  AND  GIVE  HESS 
C 

200  IC0M(ISTR+2)=IC0H(ISTR+2)+1 

IF(IC0M(ISTR+2) . LT„ IMAX)  RETURN 
URITE(LUN,205)IC0H(ISTR*2) 

ICOM(  ISTR  +  1 )  =  ICOMUSTR+5) 

RETURN 

C 

c 

220  URITE(LUN,222) 

GOTO  600 
C 

c 

300  URITE(LUN,31 0) 

GOTO  600 
C 

C*  IF  I  EQUAL  TO  ILEN,  SEARCH  FOR  UATER  COMPOSITION  FAILED. 

C 

500  CALL  EXEC( 1 1 , ITIME) 

IT=ITIME(4)*100+ITIME(3) 

URITE(LUN,  1 04) IT 

104  FORMAT ( "CANNOT  FIND  BOTTOM  COMPOSITOM  OF  UATER  AT", 15) 

C 

c 

600  IBEG  =  IADD  +  1 

IEND  =  IADD  +  ILEN 

URITE(LUN,61 2)  (ICOM(I)  ,  I  =  IBEG, IEND) 

GO  TO  200 
C 

C  FORMAT  STATEMENTS. 

C 

205  FORMAT ( '  READING  FROM  GC  HAS  NOT  CHANGE  FOR" , 15,/, 

*  '  COMPOSITION  SET  TO  DEFAULT. ' ) 

222  FORMAT ( '  GC  OUTPUT  IS  NEGATIVE'' ) 

310  FORMAT ( *  GC  READING  CHANGE  GREATER  THAN  MAX  ') 

612  F0RMAT(40A2,/,40A2) 

END 
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FTN4 

0002 

PROGRAM  GCSUT ( 3 9 80 ) ,  GC  AUTO.  SAHPLING  INITIATION  PROGRAH  D 

0003 

CALL  U A I T (5,2, IERR) 

0004 

CALL  D0L(1,1,0,4,IERR) 

0005 

IF ( IERR. ME. 1 )  URITE( 1,100)  IERR 

0006 

STOP 

0007 

100 

FORMAT ( '  GCSUT  CALL  DOL  ERROR,  IERR  =  ',12) 

0008 

END 
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FTN4 

PROGRAM  BDC99(3, 85) ,  DISTILLATION  COLUMN  MONITOR  PROGRAM  HK 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


0005 

C 

THIS  PROGRAM  MONITORS  THE 

DISTILLATION 

COLUMN  IN 

THE 

0006 

C 

ABSENCE  OF  AN  OPERATOR. 

THE  KEY  VARIABLES 

ARE: 

0007 

C 

0008 

C 

KEY  VARRIABLES  ICHAN  NO. 

I.D. 

HILIM 

LOLIM 

IDUN2 

0009 

C 

FEED  FLOU  LOOP 

12 

1 

4900 

1000 

6MIN 

0010 

C 

DISTILLATE  FLOU  LOOP 

10 

2 

4900 

1000 

30MIN 

0011 

C 

BOTTOM  FLOU  LOOP 

8 

3 

5100 

1000 

30MIN 

0012 

C 

REFLUX  FLOU  LOOP 

13 

4 

5000 

1000 

30MIN 

0013 

C 

STEAM  FLOU  LOOP 

9 

5 

3800 

1200 

15MIN 

0014 

C 

COOLING  UATER  LOOP 

11 

6 

5100 

1200 

15MIN 

0015 

C 

TOUER  PRESSURE  LOOP 

16 

7 

4900 

1200 

15MIN 

0016 

C 

REBOILER  LEVEL  LOOP 

17 

8 

4900 

1100 

15MIN 

0017 

C 

CONDENSER  LEVEL  LOOP 

18 

9 

4900 

1100 

30MIN 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

COMMON  I  COM  < 1 ) 

REAL  DCD AT  A  < 1 1 ) ,CONA< 1 1 ) , CONB<  1 1 ) 

INTEGER  IPARAM(5) , ITIME(5) , I  DATE (15) ,  ITY1 ( 1 1 ) , IT Y2 ( 1 1 ) , 
$ICHAN(11),IDATA(11),IDUN1(9),IDUN2(9),HILIM(9),L0LIM(9) 

DATA  ICHAN/1 2, 10, 8, 13, 9, 11, 16, 17, 18, 15,0/ 

DATA  CONA/2. 5777, 1.5205, 1.51 05, 2. 2975, 2. 25, 71. 278, 
$0.96988,1.0212,1.0723,0.12857,0.0010/ 

DATA  CQNB/9*0. 0,89. 167,0.0/ 

DATA  NDATA,NTEST/1 1 ,9/,ITY1/10*1 ,2/, ITY2/8*1 ,3*2/ 

DATA  IDUN 1/9*0/, IDUN2/5,3*30,4*1 5,30/,IDQUN/0/ 

DATA  HILIM/2*4900,51 00,3800,5000,51 00,3*4900/ 

DATA  L0LIM/4*1 000, 3*1200, 2*1100/ 


C 

C 

C 

2000 

2010 


2020 

2030 

C 

C 

C 


FORMAT  STATEMENTS 

FORMAT </,5X,**BDC99  UP  8  ”,15A2) 
F0RMAT(/,1X,"TIME“,2X,"C.L.  ERR", IX,' 
$"  3  -  BP  ",1X,H  4  -  RE  ",1X,"  5  -  ST 
$H  7  -  PR  H,1X,M  8  -  RL  W,1X,"  9  -  CL 
$-11  -  XB  M) 

FORMAT ( 1Xv2I2y12( 1X,F8.3) ) 

FORMAT ( /// ) 

INITIALIZATION 

CALL  RMPAR(IPARAM) 

LUN=IPARAM( 1 ) 

MULT=IPARAH(2) 

I RES-3 

I CHAN ( 1 1 >  * ICOM ( 8 ) ♦ 1 
DO  5000  1=1,9- 
IDUN2(I)=IDUN2< I) /MULT 


1  -  FE  '• 
",1X,"  6 
”  1  X , H 1  0  - 


,1X,' 

-  cu 

XD  ' 


2  - 
",1X, 
,1X, 


TP 


■ 


. 


0052 

0053 

0054 

0055 

0056 

0057 

0058 

0059 

0060 

0061 

0062 

0063 

0064 

0065 

0066 

0067 

0068 

0069 

0070 

0071 

0072 

0073 

0074 

0075 

0076 

0077 

0078 

0079 

0080 

0081 

0082 

0083 

0084 

0085 

0086 

0087 

0088 

0089 

0090 

0091 

0092 

0093 

0094 

0095 

0096 

0097 

0098 

0099 

0100 

0101 

0102 


203 


5000  CONTINUE 

CALL  FT  I HE ( I DATE ) 

URITEt 1 ,2000)  IDATE 
URITE(LUNr2000)  IDATE 
URITE(LUN,2010) 

I  SKI P—0 
C 

C  DATA  ACQUISITION 
C 

5010  CALL  ACQUI<LUN,NDATA,ITY1,ITY2,ICHAN,IDATA,DCDATA,C0NA,C0NB, 
♦IERR) 

IF  ( IERR.NE. 1 )  GO  TO  5005 
CALL  EXEC ( 1 1 , IT  I ME ) 

ERR=100.0*(DCDATA( 1 )-DCDATA(2)-DCDATA(3) )/BCDATA(1 ) 
URITE(LUN,2020)  ITIME(4) ,ITIME(3) ,ERR, (DCDATA( I) ,  1=1 ,NDATA) 
C 

C  VARIABLE  LIMIT  CHECKING 

C 

CALL  LMTCK(LUN,NTEST, IDATA,HILIM,LOLIH, IDUN1 , IDUN2, IDOUN, 
SISKIP) 

IF ( IDOUN.NE.O)  GO  TO  5005 
C 

C  SUSPEND  THE  PROGRAM  FOR  MULT-MIN. 

C 

CALL  UAIT(MULT, IRES, IERR) 

IF  ( IERR.NE.  1 )  GO  TO  5005 
ISKIP=ISKIP+1 

IF  ( ISKIP.NE.45)  GO  TO  5010 
ISKIP=0 

URITE(LUN,2030) 

URITE(LUN,201 0) 

GO  TO  5010 
C 

C  DISTILLATION  COLUMN  SHUT  DOUN 
C 

5005  CONTINUE 

CALL  DOL( 1 , 1 ,8,8, IERR) 

C 

C  RESET  THE  ECO  15  SEC.  AFTER  THE  SHUT  DOUN 
C 

CALL  WAIT ( 15, 2, IERR) 

CALL  DOL( 1 , 1 ,0,8, IERR) 

STOP 

END 

C 

C 

SUBROUTINE  LMTCK(LUN,NDATA,IDATA,HILIM,LOLIM, IDUN1 ,IDUN2,ID0 
SISKIP) 

INTEGER  IDATA(1),HILIM(1),L0LIM(1),IDUN1 <1 ), IDUN2( 1 ) , ITIME( 5 
♦IDATE(15),LABEL(9) 

DATA  LABEL/2HFE,2HTP,2HBP,2HRE,2HST ,2HCU,2HPR,2HRL,2HCL/ 


. 


0103 

2000 

FORMAT !/,5X,HA  L  E  R  T  LOOP  -  -,A2) 

0104 

2010 

FORMAT !/,5X, HD  ANGER  DIST •  COL.  DOUN  0  " 

0105 

$MDOUN  LOOP  =",A2) 

0106 

2020 

FORMAT !/,5X,"BDC99  DOUN  8  ", 15A2,5X, "D.L.  =  M,A2 

0107 

CALL  FTIME(IDATE) 

0108 

DO  5000  1=1 r NDATA 

0109 

IF (IDATA(I) .LE.LOLIM< I ) .OR. 

0110 

$IDATA(I) .GE.H I L I M < I ) )  GO  TO  5010 

0111 

GO  TO  5020 

0112 

5010 

IF ( IDUN1 ( I ) .GT. IDUN2! I) )  GO  TO  5030 

0113 

IDUN1 < I )=IDUN1 ( I )+1 

0114 

ID0UN=0 

0115 

URITE(LUN,2000)  LABEL ( I ) 

0116 

ISKIP=ISKIP+3 

0117 

GO  TO  5000 

0118 

5020 

IDUN1 < I >  =0 

0119 

ID0UN=0 

0120 

5000 

CONTINUE 

0121 

RETURN 

0122 

5030 

I D  OUN= 1 

0123 

URITE!LUN,201 0)  IDATE,LABEL! I) 

0124 

WRITE! 1  ,2020)  I DATE p LABEL ( I ) 

0125 

RETURN 

0126 

END 

0127 

C 

0128 

C 

0129 

SUBROUTINE  ACQUI!LUN,NDATA,ITY1 , IT Y2, I CHAN, IDATA 

0130 

$CONA,CONB,IERR) 

0131 

COMMON  I COM ( 1 ) 

0132 

DIMENSION  DCDATA(1),C0NA(1),C0NB!1) 

0133 

INTEGER  ICHANd ) , I DATA ( 1 ) , ITY1 ( 1 ) ,ITY2! 1 ) 

0134 

2000 

FORMAT!/, 5X,MDIST.  COL.  -  LOOP  ",I2,"  IERR 

0135 

DO  5000  1=1, NDATA 

0136 

IF  ! ITY1 ! I ) .EQ.2)  GO  TO  5010 

0137 

CALL  AIRDI1, ICHANd), ID  AT  A(I), IERR) 

0138 

IF  ! IERR.EQ. 1 )  GO  TO  5020 

0139 

URITE(LUN,2000)  I, IERR 

0140 

RETURN 

0141 

5020 

IF  ! IDATA! I) . GE - 1 000 . AND . I DATA < I ) .LE.5000)  GO  TO 

0142 

IF  ! IDATA! I ) .LT. 1 000)  IDATA! I ) =1 000 

0143 

IF  ! IDATA! I ) .GT.5000)  IDATA! I ) =5000 

0144 

5021 

RDATA=0.0250*FL0AT( IDATA ! I ) ) -25.0 

0145 

GO  TO  5030 

0146 

5010 

IERR=1 

0147 

IDATA! I )  =  ICOM! ICHAN! I ) ) 

0148 

RDATA=FLOAT ( IDATA! I ) ) 

0149 

GO  TO  5030 

0150 

5030 

IF  ( ITY2! I ) .EQ.2)  GO  TO  5040 

0151 

DC DAT A! I )=CONA< I )*SQRT  <RDATA)+CONB! I ) 

0152 

GO  TO  5000 

0153 

5040 

DCDATA! I )=CONA( I )*RDATA+CONB! I ) 

0154 

5000 

CONTINUE 

0155 

RETURN 

0156 

END 

■ 


